% Gauss version written by: Gautam Gowrisankaran, April 25, 1993
% Modified by: Ron Goettler, Feb. 27, 2013
% converted to Matlab mid 2013.
% This program generates the welfare output programs used for the Markov-perfect game.
function simstuff = welf_ma(batchnum)
disp('entering welf_ma.m ...');
disp('**** Computing Welfare Statistics ****');

% Build Par struct with general program data
Par.batchjob        = batchnum;
Par.configfile      = ['config_' num2str(Par.batchjob) '.mat'];
Par.small           = 1e-8;
Par.nocheck         = -1e50; 
Par.uu_competition  = 1;
Par.uu_monopoly     = 2;
Par.uu_planner      = 3;
Par.uu_quality      = 1;
Par.uu_cost         = 2;
Par.uu_capacity     = 3;
Par.uu_ran_entry    = 1;
Par.uu_det_entry    = 2;

% Load Params struct from file (if it exists) or abort
if exist(Par.configfile, 'file')
    disp('loading config file...');
    load(Par.configfile, 'Params');
else
    disp(['Could not find config file ' Par.configfile]);
    disp('Aborting!');
    return;
end

% Display model type
[prefix,eqltype,invtype] = acronym(Par, Params.eql_type, Params.ind_type);
disp(['model: ' eqltype ', investment in ' invtype]);

% Define local variables
M  = Params.mkt_size;
mc = Params.mc;
beta = Params.beta;
ParRlg.kmax=Params.kmax;
ParRlg.no_force_exit=0;
if ((Params.entry_high>1e99 && Params.scrap_val<-1e99) || (Params.entry_high<.0001 && Params.scrap_val<.0001))
    ParRlg.no_force_exit = 1;
end


% Set the default random number generator stream
RandStream.setGlobalStream( RandStream.create('mrg32k3a','Seed',58921) );

wstart = zeros(Params.max_firms,1);

if (Params.entry_at > 0)
    entry_T = Params.entry_at;                  % entry at fixed rung on quality ladder
else
    entry_T = ParRlg.kmax+Params.entry_at;      % entry at  abs(Params.entry_at)  below  leader
end

if (Par.batchjob==0)
    disp('Enter initial efficiency levels.  Defaults are zero (meaning out).');
    disp(['Valid range for parameter values is 0 - ' int2str(ParRlg.kmax)]);
    for i=1:Params.max_firms
      wdef=0;
      if i==1; wdef= entry_T+2;  end;
      prompt=['Firm ' int2str(i) ': '];
      wstart(i)=getint(Par, prompt, 0, ParRlg.kmax, wdef);
    end
    numtimes=getint(Par, 'Number of periods to simulate [1-1K, default 100]: ',1,1000,100);
    numruns =getint(Par, 'Number of runs to simulate for [1-100K, default 1000]: ',1,100000,1000);
else
    wstart(1) = ParRlg.kmax;                                               % AFTER load isentry[], wstart(2) = entry_T if isentry>0 with this monop
    if (Params.wstar>0)
        wstart(1) = min([ParRlg.kmax entry_T+2]);                          % alt: minc(kmax-1|(ceil((6+pp[_WSTAR]-Params.rlg_wshift)/Params.rlg_wscale)))  where WSTAR is on post shift/scale grid, but wstart is on 0,1,2...Kmax grid
    end

    if (Params.max_firms==2 || Params.rlg_y==0)
        wstart(2) = entry_T;                                                % fixed DUOP, so start with DUOP, or linear model(since want to avoid monopoly)
    end
    if ((Params.entry_high<1e-8 || Params.entry_high>1e100) && Params.wstar<=0)
        if 0
            wstart = ParRlg.kmax+zeros(Params.max_firms,1);                % --> start with all firms at frontier,  or
        else
            wstart = [ParRlg.kmax ; (entry_T+zeros(Params.max_firms-1,1))];% --> start with all laggards at entry_T
        end
    end
    numtimes = 100;   
    if (Params.wstar>0)
        numtimes= 500;                                                      % increase numtimes if WSTAR to remove uncertainty of whether wstart(1) is near wfinal(1).  If WSTAR<=0  wstart(1) = wfinal(1) = kmax
    end
    numruns = 1000;
end

% Build binomial table
binom = eye(Params.max_firms+ParRlg.kmax+1);
binom = [zeros(Params.max_firms+ParRlg.kmax+1, 1) binom];
for i = 2:Params.max_firms+ParRlg.kmax+1;
    binom(i, 2:i) = binom(i-1, 2:i)+binom(i-1, 1:i-1);
end
if (Params.wstar == 0)
    ParRlg.kmax = ParRlg.kmax-1; 
end

wmax = binom(Params.max_firms+ParRlg.kmax+1,ParRlg.kmax+2);       % Number of possible industry structures


% Load in all the data stored by the equilibrium generation program 
% This data is: v (value), x (investment), p (probability of state rising),
% isentry
filename = [prefix 'markov.' int2str(Params.max_firms) 'ot' '.mat'];
load(filename);
v = newvalue;           % value functions
x = newx;               % investment policy
p = probw;              % prob(own w up)
isexit = newexit;       % prob exit for random scrap
%isentry = isentry;      % prob entry
v_entry = oldv_entry;   % value of entry


% Load in all data from the static profit calculation
% The data is: firm profits, consumer surplus, market shares at each state,
% price/cost margins, one-firm concentration ratios, prices.
filename = [prefix 'cons.' int2str(Params.max_firms) 'f' '.mat'];
load(filename);
% RP - profit is agprof here, consistent with welf_ma.g, but is it correct?
profit = agprof;        % firm profits
%csurplus = csurplus;    % consumer surplus
%share = share;          % market shares at each state
pmcmargm = pmcmarg;     % price/cost margin
concentm = concent;     % concentration ratios
%price = price;          % prices


if (Par.batchjob>0 && isentry(encode(binom,Params,wstart))>0)
    wstart(2) = entry_T;
end
wstart = sortrows(wstart,-1);

disp('-----  Simulation of Dynamic Oligopoly  ---------------------------------------------------------');
disp(['Starting conditions: [' num2str(wstart') ']']);
disp(['Number of simulation periods: ' int2str(numtimes)]);
disp(['Number of runs: ' int2str(numruns)]);

%if (Params.wstar==0)
%    ParRlg.kmax = ParRlg.kmax-round(2+Params.kmax/6*(1-Params.delta));
%    disp(['** Adjusting kmax down by ' round(2+Params.kmax/6*(1-Params.delta)) ' (to ' ParRlg.kmax ') in welf_ma since using WSTAR = 0 is more accurate when simulate with kmax lower than used to get investment policy']);
%end

wmax = binom(Params.max_firms+ParRlg.kmax+1,ParRlg.kmax+2);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Define matrices to hold all simulation results
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
consurp  = zeros(numruns,1);    % (Total) consumer surplus
prodsurp = zeros(numruns,1);    % Producer surplus

prodsurp2= zeros(numruns,1);    % PS with unbounded outside utility using reoptimized prices  
consurp2 = zeros(numruns,1);    % CS with unbounded outside utility using reoptimized prices  

totinnov = zeros(numruns,1);    % other RLG stats   
totinnov1= zeros(numruns,1);    % leader innovation 
totinnov2= zeros(numruns,1);    % firm 2 innovation 
totinnovT= zeros(numruns,1);    % frontier innov T  
totinvest= zeros(numruns,1);    % industry invest   
totinves1= zeros(numruns,1);    % leader   invest   
totinves0= zeros(numruns,1);    % industry invest 0 && no entry or exit --> Absorbing State 
totinve0n= zeros(numruns,1);    % n firms invest 0  
totprof1 = zeros(numruns,1);    % leader profit excluding investment 
totprof2 = zeros(numruns,1);    % firm 2 profit excluding investment 
totnfirms= zeros(numruns,1);
tot_firms= zeros(numruns,Params.max_firms+1);  % last col is freq of empty industry
totentry = zeros(numruns,1);
totshare0= zeros(numruns,1);
totshare1= zeros(numruns,1);
totshare2= zeros(numruns,1);
totprice1= zeros(numruns,1);
totprice2= zeros(numruns,1);
totw2    = zeros(numruns,1);
totw2T   = zeros(numruns,1);
totshar1T= zeros(numruns,1);
totmarkup= zeros(numruns,1);
totwstar = zeros(numruns,1);
totwsame = zeros(numruns,1);
totLRstay= zeros(numruns,1);    % a firm at Lowest Rung stays with some probability
totLRx   = zeros(numruns,1);    % a firm at Lowest Rung invests > 0 
totpdiff = zeros(numruns,1);    % static p + RLGforced  -  newton p 

totexit  = zeros(numruns,1);    % total number of firms that exit, and therefore contribute to totv0_vx 
totv0    = zeros(numruns,1);    % initial v 
totvx    = zeros(numruns,1);    % realized profits 
totv0i   = zeros(numruns,1);    % initial v, initial leading firm only 
totvxi   = zeros(numruns,1);    % realized profits, initial leading firm only 
totlspi  = zeros(numruns,1);    % lifespan, initial leading firm only 
totxT0   = zeros(numruns,1);    % # runs ending in a state with zero investment (absorbing state if delta=0) 
discEU_T = zeros(numruns,1);    % discounted EU from final period (check if small as needed for run to be long enough and CS to be well-defined) 
discCS_T = zeros(numruns,1);    % discounted CS from final period (check if small as needed for run to be long enough and CS to be well-defined) 
hitkmax  = zeros(numruns,1);    % # times lead firm at kmax 
maxk     = 0;                   % maximum k achieved across all runs 
x1absorb = 0;                   % leader investment in absorbing state in which leader invests to deter laggards' investments 

active   = zeros(Params.max_firms+1,1);  % No. of periods with n firms active.
exitors  = 0;  % No. of exitors
entrants = 0;  % No. of entrants
entexit  = 0;  % No. of periods with both entry and exit 
lifedis  = 0;  % Distribution of firm lifespans 
valuedis = 0;  % Distribution of firm's total profits 


%{
simfile = "BJ" $+ ftos(BatchJob,"%*.*lf",1,0) $+ "_sim";
output file = ^simfile reset;   screen off;    output on;
"% nr(1)  t(2)  cum_innov(3)  pr_innov(4)  pr_entry(5)   w  thisexit   ax/(1+ax)   share  price  profit";
screen on;     output off;
%}


for nr=1:numruns
    RLGforced = 0;
    RLGinnov  = 0;
    discEU    = zeros(numtimes,1);      % discounted EU : used to compute CS when RLG_out = 0 and Params.rlg_y > 0
    discCS    = zeros(numtimes,1);      % discounted CS : used to compute CS when RLG_out = 0 and Params.rlg_y > 0

    firm_id  = 1;
    if Params.max_firms>1;
        firm_id  = [1 ; zeros(Params.max_firms-1,1)];       % track initial value and realized values of 1st leading firm
    end
  
    lifemx   = zeros(Params.max_firms,1);  % No. of periods that firm has been active for
    valuemx  = zeros(Params.max_firms,1);  % Total profits of each firm to date
    valuem0  = zeros(Params.max_firms,1);  % value of each firm's initial state
    wthis    = wstart;
    wthisprev= wstart;
    lifemx   = lifemx + (wthis > 0);

    ii = (lifemx==1);  % initialize with firms that already exist prior to 1st simulated period
    if (sum(ii)>0)
        valuem0 = valuem0.*(1-ii) + ii.*(v(encode(binom,Params,wthis),:)');
    end
    
    for t=0:(numtimes-1)
        codew2 = encode(binom,Params,wthis);    % Get probabilities of investment causing a rise in eff level & actual investment and value function for this state
        prob = p(codew2,:)';                    % prob(own w up)
        xx = x(codew2,:)';                      % Investment decisions made prior to observing entry/exit outcomes
        vv = v(codew2,:)';                      % Profit should also be based on codew2 since also before outcomes

        if ((t+1==numtimes) && mod(nr,(numruns/10)) == 0)   % print last state for 10 of the runs -- shows progress in completing numtimes, and gives limited sample of distribution of terminal state
            disp(['Runs computed: ' sprintf('%6.0f',nr) '  w_T = ' sprintf('%2.0f',wthis')  '  invest = ' sprintf('%8.3f',xx') '  innov = ' sprintf('%7.4f',prob') '  entry = ' sprintf('%7.4f',isentry(codew2)) '  lifemx = ' sprintf('%6.0f',lifemx') '  valuemx = ' sprintf('%8.1f',valuemx') '  valuem0 = ' sprintf('%8.1f',valuem0') '  Outside Good: #Innovation, #Free Bumps = ' sprintf('%6.0f',[RLGinnov RLGforced]) '  discEU avg~final = ' sprintf('%12.1f', [mean(discEU) discEU(numtimes)]) ]);
        end

        if (sum(wthis==1 & xx>1e-8)>0)
            totLRx(nr) = totLRx(nr) + 1;   
        end
        if (sum(wthis==1 & isexit(codew2,:)'<(1-1e-8))>0)
            totLRstay(nr) = totLRstay(nr) + 1;
        end
    
        ii = sum(wthis>0);
        if (ii==0)
            tot_firms(nr,Params.max_firms+1) = tot_firms(nr,Params.max_firms+1) + 1;
            if (sum(tot_firms(1:nr,Params.max_firms+1)) < 5)
                disp(['HEY: no firms are in the industry !!   t = ' num2str(t) ' ; wthis = ' num2str(wthis') ' ; wthisprev = ' num2str(wthisprev') ' at row ' num2str(encode(binom,Params,wthisprev)) ' in policy functions']);
            end
        else
            tot_firms(nr,ii)= tot_firms(nr,ii) + 1;                 % based on pre- entry/exit/investment   
            totnfirms(nr) = totnfirms(nr) + ii;                     % before exit & entry since profit & CS 
            totinve0n(nr) = totinve0n(nr) + sum(xx(1:ii)==0);       % pre-exit to track deadbeat firms. since totinve0n focuses on laggards, ignore RLG_minx1 as the "zero" for frontier firms 
        end
    
        totinves1(nr)= totinves1(nr) + xx(1);
        totinnov1(nr)= totinnov1(nr) + prob(1);
        totinnov2(nr)= totinnov2(nr) + prob(2);
        totprof1(nr) = totprof1(nr)  + profit(codew2,1);
        totprof2(nr) = totprof2(nr)  + profit(codew2,2);
    
        totw2(nr) = totw2(nr) + wthis(1)-wthis(2);                  % leader advantage
        
        maxk = max([maxk ; wthis(1)]);
        if (wthis(1)==ParRlg.kmax)
            hitkmax(nr)  = hitkmax(nr)+1;
        end
        if (Params.wstar>0)
            totwstar(nr) = totwstar(nr) + sum((wthis.*Params.rlg_wscale+Params.rlg_wshift)>Params.wstar);
        end
        if (Params.max_firms>1)
            totwsame(nr) = totwsame(nr) + (wthis(1)==wthis(2));
        end
        
        % Find out which firms want to leave.
        wtrans = zeros(Params.max_firms,1);
        if (Params.rlg_scrap_hi==0)                                         % fixed scrap
            [min_vv, mind_vv] = min(vv);
            ii = (min_vv == Params.scrap_val)*(mind_vv-1) + (min_vv > Params.scrap_val)*Params.max_firms;
            if (ii > 0)
                wtrans(1:ii) = wthis(1:ii);                                 % ii = #firms remaining
            end    
        else                                                                % random scrap
            prexit = rand(Params.max_firms,1);
            wtrans = wthis .* (prexit > isexit(codew2,:)');                 % exit if prexit <= isexit
        end

        % Now figure out exit. Must capture firms who exit voluntarily,
        % as well as firms whose efficiency has been driven down to zero.
        lifemx = lifemx + (wthis > 0);
        thisexit = (wtrans == 0) & (lifemx > 0);                            % checking lifemx ensures the exitor is a real firm, not a vacancy
        lifemx = lifemx - (wthis > 0);

        phis = zeros(Params.max_firms,1);                                   % vector of realized Params.scrap_val values.  Want simulated phis, not E(phi|exit)
        if (sum(thisexit) > 0)
            ii = find(thisexit==1);
            if (Params.rlg_scrap_hi==0)
                phis(ii) = Params.scrap_val*ones(size(ii));                 % ii = exitor indices
            else
                phis(ii) = Params.rlg_scrap_hi - (Params.rlg_scrap_hi-Params.scrap_val)*prexit(ii);  % prexit are the rand() draws compared to isexit() above
            end
            
            xx   =   xx .* (1-thisexit);    % zero out exitors xx and prob since those values are conditional on NOT exiting
            prob = prob .* (1-thisexit);    % only relevant for random scrap case, but harmless to do for all cases
        end
        
        % codew = encode(binom,Params,wtrans); HEY: why does PM original code use wtrans (i.e., AFTER EXIT market structure)??
        codew = codew2;             % eql_ma_foc.m uses profit(w,j) where w is BEFORE exit, so RLG uses codew2 NOT encode(Par,Params,wtrans)
        sigma = share(codew,:)';    % may get overwritten later

        
        % Now, tally the statistics.  RLG moved this block from post-entry
        % determination to here (since profits are BEFORE entry/exit/invest
        % outcomes
        consurp(nr) = consurp(nr) + (beta^t) * csurplus(codew);                                                    % if Params.rlg_y > 0 will convert to true CS after last period
        consurp2(nr)= consurp2(nr)+ (beta^t) *(csurplus(codew)+M*RLGforced*Params.rlg_wscale/Params.rlg_alpha);    % add extra CS (or M/alpha scaled EU) if grid has shifted
        discEU(t+1) =               (beta^t) *(csurplus(codew)*Params.rlg_alpha/M + RLGforced*Params.rlg_wscale);  % enables checking whether EU and CS from last period is small after being discounted to 0

        if (Params.rlg_y)
            discCS(t+1) = (beta^t) * M * ( exp( (discEU(t+1)/(beta^t))/Params.rlg_alpha ) - Params.rlg_y );
        else
            discCS(t+1) = (beta^t) *(csurplus(codew)+M*RLGforced*Params.rlg_wscale/Params.rlg_alpha);
        end

        i = sum(wthis>0);       % #active firms
        if (i==0)               % EMPTY industry
            prodsurp(nr) = prodsurp(nr)  + (beta^t)*(sum(phis) - sum(xx));
            prodsurp2(nr)= prodsurp2(nr) + (beta^t)*(sum(phis) - sum(xx));
            valuemx = valuemx + (beta.^(lifemx-1)).*(-xx);                   % NOTE: firm that exists in wstart has lifemx = 1 when t=0, so use lifemx-1
      
            totinves0(nr)= totinves0(nr) + (isentry(codew)==0);             % absorbing state
            totinvest(nr)= totinvest(nr) + sum(xx);                         % consurp2, totmarkup, totshare0, are unchanged
            totshare0(nr)= totshare0(nr) + 1;
        else                    % industry NOT EMPTY
            sigma = share(codew,1:i)';
            ppp   = price(codew,1:i)';

            www = M * (ppp-mc).*sigma - Params.fixed_cost;              % current period profit
            prodsurp2(nr)= prodsurp2(nr) + (beta^t)*(sum( www )            + sum(phis) - sum(xx)) ; % entryfee subtracted later
            prodsurp(nr) = prodsurp(nr)  + (beta^t)*(sum(profit(codew,:)') + sum(phis) - sum(xx)) ;

	    valuemx = valuemx + (beta.^(lifemx-1)).*( profit(codew,:)' - xx );      % NOTE: firm that exists in wstart has lifemx = 1 when t=0, so use lifemx-1.  Enables monopolist valuemx to match its prodsurp

            totinvest(nr)= totinvest(nr) + sum(xx);
            if (Params.delta==0)
                totinves0(nr)= totinves0(nr) + (sum(xx)==0 && isentry(codew)==0 && sum(isexit(codew,:)'.*wthis)==0 );   % Absorbing State
                if (ParRlg.no_force_exit && xx(1)>0 && max(wthis(2:Params.max_firms))<=1)  %  Can also have absorbing state with Sole Frontier x > 0 and all other firms at 1 or out and all entry, exit = 0, IF ParRlg.no_force_exit = 0
                    totinves0(nr)= totinves0(nr) + (sum(xx(2:Params.max_firms))==0 && isentry(codew)==0 &&  sum(isexit(codew,:)'.*wthis)==0);   % Absorbing State, but with investment
                    x1absorb = xx(1);
                end
            end
            totmarkup(nr)= totmarkup(nr) + sum(ppp.*sigma) / mc / sum(sigma);
            totshare0(nr)= totshare0(nr) + 1-sum(sigma);
            totshare1(nr)= totshare1(nr) + sigma(1);
            totprice1(nr)= totprice1(nr) + ppp(1);
	
            if (i>1)
                totshare2(nr)= totshare2(nr) + sigma(2);  
                totprice2(nr)= totprice2(nr) + ppp(2);  
            end
            %end % end of unbounded state space if, which was commented out
        end  % end of block moved by RLG to precede exit instead of follow entry stuff
        
        % Handle exiting firms
        if (sum(thisexit) > 0)
            % valval = selif(valuemx,thisexit)+(beta.^( @ 1+ @ selif(lifemx,thisexit)).*selif(phis,thisexit));
            valval = selif(valuemx,thisexit)+(beta.^(selif(lifemx,thisexit)).*selif(phis,thisexit)); % RLG added profit() here since exit occurs AFTER profits
            
            if ((firm_id' * thisexit )>0)               % the leading firm exited
                i = find(firm_id==1);
                totlspi(nr)= lifemx(i);                 % mean totv0i = mean totvxi  if value function correct
                totv0i(nr) = valuem0(i);
                totvxi(nr) = valuemx(i) + beta.^lifemx(i)*phis(i);
                firm_id = zeros(Params.max_firms,1);    % no longer have the initial leading firm
            end
      
            % lifedis  = [lifedis ; selif(lifemx,thisexit)];
            % valuedis = [valuedis ; valval];
            totexit(nr)= totexit(nr) + sum(thisexit);
            totv0(nr)  = totv0(nr) + sum(selif(valuem0,thisexit));
            totvx(nr)  = totvx(nr) + sum(valval);
            lifemx     = lifemx  .* (1-thisexit);
            valuemx    = valuemx .* (1-thisexit);
            valuem0    = valuem0 .* (1-thisexit);
        end

        
        % Handle entering firm
        entrypr = rand(1,1);           % original PM used  encode(binom,Params,wtrans), which is AFTER exit occurs, in place of codew below
        if entrypr < isentry(codew);   %{ with random scrap, entry/exit DECISIONS both depend on wthis NOT wtrans (i.e., BEFORE either entry or exit is observed)
            % prob(Params.max_firms)  = 0;
            % when nfirms > 2, only do leap if TWO open spots so can still check frequency of full industry for need to increase nfirms
            doleap = (Params.rlg_leap>0 && Params.max_firms>1);
            if (doleap && Params.max_firms>2)
                doleap = (wthis(Params.max_firms-1)==0);
            end
            if (doleap && rand(1,1) < Params.rlg_leap)
                wtrans(Params.max_firms) = ParRlg.kmax;
            else
                if (Params.entry_at>0)
                    wtrans(Params.max_firms) = Params.entry_at;
                else    % enter at |Params.entry_at| below start-of-period frontier firm even if not at kmax.  enter relative to kmax if empty industry
                    wtrans(Params.max_firms) = max([1;(wthis(1)+Params.entry_at)]) + (wthis(1)==0)*(ParRlg.kmax+Params.entry_at-1);
                end
            end
            entryfee     = Params.entry_low + entrypr * (Params.entry_high - Params.entry_low);   % uses realized entryfee, NOT E(fee|enter)
            totentry(nr) = totentry(nr)+1;
            valuemx(Params.max_firms) = 0;
            valuem0(Params.max_firms) = v_entry(codew);

            prodsurp(nr)  = prodsurp(nr)  - (beta^t)*entryfee;    % PM subtracted entryfee earlier
            prodsurp2(nr) = prodsurp2(nr) - (beta^t)*entryfee;    % RLG moved here so that profit earned before entry/exit/invest
        end % entry
    
        lifemx = lifemx + (wtrans > 0);  % moved here from before profit tallied by RLG since want lifemx = 0 for new entrant

        
        % RLG moved wnext stuff AFTER computing surpluses since do not want to update RLGforced until after computed CS
        wnext = wtrans+(prob>=rand(Params.max_firms,1));
        if (rand(1,1) <= Params.delta)
            if ParRlg.no_force_exit;
                wnext = wnext-1 + (wnext==1);  % Never bump firms off lowest rung of ladder
            else
                wnext = wnext-1;
            end
            RLGinnov = RLGinnov + 1;
        end
        if (Params.wstar <= 0 && sum(wnext==ParRlg.kmax+1)>0)  % RLG statespace && frontier moving beyond kmax, so shift ALL down
            if (ParRlg.no_force_exit)
                wnext = wnext-1 + (wnext==1);  % Never bump firms off lowest rung of ladder
            else
                wnext = wnext-1;
            end
            RLGforced = RLGforced+1;
        end
        wnext = max([wnext zeros(Params.max_firms,1)]')';
        wnext = min([wnext ParRlg.kmax*ones(Params.max_firms,1)]')';        % RLG added

        %if nr<101;  %{ save to file for plotting by  simplots.m :  nr  t  cum_innov  pr_innov  pr_entry  w  thisexit  ax/1+ax  share  price  profit %}
        %    output file = ^simfile;   screen off;    output on;           %{ pr_innov uses investment ignoring exit,   fine since leader rarely exits %}
        %    format /m1 /rds 1,0;  nr  t+1  RLGforced+RLGinnov+wthis(1)-wstart(1);;
        %    format /m1 /rds 6,4;  (1-(1-p(codew,1))^sum(wthis.==wthis(1)))  isentry(codew) ;;
        %    format /m1 /rds 1,0;  wthis' thisexit';;
        %    format /m1 /rds 6,4;  p(codew,.)  share(codew,.)  price(codew,.)  profit(codew,.) ;
        %    screen on;     output off;
        %end ; 

        % Now re-sort all firm level data, to reflect the fact that firms must be in descending order next year
        wthisprev = wthis;
        temp = sortrows([wnext lifemx valuemx valuem0 firm_id],-1);
        wthis   = temp(:,1);  
        lifemx  = temp(:,2);  
        valuemx = temp(:,3);  
        valuem0 = temp(:,4);  
        firm_id = temp(:,5);
    end % numtimes repetitions

    % get realized + cont.value of firms active in last simulated period
    thisexit = (lifemx>=1);  % RLG confirmed valuemx = valuem0 exact when wstart is an absorbing state (so discounting powers correct)
    if (sum(thisexit)>0)     % RLG also notes that mean totvx usually > mean totv0 ... probably due to value driven by outliers?
        valval = selif(valuemx,thisexit) + beta.^(selif(lifemx-1,thisexit)).*selif(v(encode(binom,Params,wthis),:)', thisexit);
        totexit(nr) = totexit(nr)  + sum(thisexit);
        totv0(nr) = totv0(nr) + sum( selif(valuem0,thisexit)  );  % mean totv0 = mean totvx  if value function correct
        totvx(nr) = totvx(nr) + sum( valval );
    end
    
    if (sum(firm_id==1)>0)
        totv0i(nr) = selif(valuem0,firm_id==1);
        totvxi(nr) = selif(valuemx,firm_id==1) +  beta^numtimes * selif( v(encode(binom,Params,wthis),:)', firm_id==1);
        totlspi(nr)= selif( lifemx,firm_id==1);  % will be numtimes+1
    end
    totxT0(nr) = (totinves0(nr)>0);  % any periods hit absorbing state
    totinnov(nr) = RLGforced+RLGinnov+wthis(1)-wstart(1);

    % Conceptually 2 ways to compute dynamic CS when utility has log(y-p): 
    % see discussion in profit_main.m. 
    % discCS(t+1)= (beta^t) * M * exp( ln( sum( exp( (RLGforced + wthis(ii)).*Params.rlg_wscale+Params.rlg_wshift + Params.rlg_alpha * ln( Params.rlg_y - price(codew,ii)' ) )))/Params.rlg_alpha);
    % alpha*ln(y)/(1-beta) = consurp2 = cumsum( discounted EU ) --> y as a per-period compensating income = exp( (1-beta)*consurp2/alpha ) --> Y= y/(1-beta) is the discounted compensating y.  
    % below is only valid if periods per run is large, since using 1/(1-beta) which is infinite sum 
    % meanc(discEU())*beta^numtimes/(1-beta)) assumes all periods BEYOND numtimes have this EU, which is still too low, but better than assuming EU = 0 after numtimes
    if (Params.rlg_y > 0)
        www = (1-beta)*consurp(nr)*Params.rlg_alpha/M;                      % find stream of income that generates this per-period equivalent EU to lifetime discounted EU
        consurp(nr) = M * ( exp( www/Params.rlg_alpha ) - Params.rlg_y) /(1-beta);
  
        www = (1-beta)*consurp2(nr)*Params.rlg_alpha/M;
        consurp2(nr)= M * ( exp( www/Params.rlg_alpha ) - Params.rlg_y) /(1-beta);
    end
  
    discEU_T(nr) = discEU(numtimes);
    discCS_T(nr) = discCS(numtimes);

    if (sum(consurp<0)>0)
        disp('HEY: in welf_ma.m CS <0,      resetting to 0');  
        consurp = max([zeros(1,numruns) ; consurp'])';    
    end
    if (sum(consurp2<0)>0)
        disp('HEY: in welf_ma.g CS2<0,      resetting to 0');
        consurp2= max([zeros(1,numruns) ; consurp2'])';   
    end
    if (sum(discCS_T<0)>0)
        disp('HEY: in welf_ma.g discCS_T<0, resetting to 0');
        discCS_T= max([zeros(1,numruns) ; discCS_T'])';
    end
    
    totw2T(nr)    = wthisprev(1)-wthisprev(2);                      % leader advantage in last period
    totshar1T(nr) = sigma(1);                                       % leader share in last period
    totinnovT(nr) = 1-(1-prob(1))^sum(wthisprev==wthisprev(1));     % prob frontier innovates last period.  Assumes no firms at frontier will exit.
end % numruns repetitions


disp(' ');
disp(['Total Enter ~ Exit (includes firms alive at end) ~ Exit-Entry-initfirms = ' num2str(sum(totentry)) ' ~ ' num2str(sum(totexit)) ' ~ ' num2str(sum(totexit-totentry)-numruns*sum(wstart>0))]);
disp(' ');

totinnov = totinnov /numtimes ;
totinnov1= totinnov1/numtimes ;
totinnov2= totinnov2/numtimes ;  % uses 0 for states in which only 1 firm
totinvest= totinvest/numtimes ;
totinves1= totinves1/numtimes ;
totinves0= totinves0/numtimes ;
totinve0n= totinve0n/numtimes ;
totnfirms= totnfirms/numtimes ;
tot_firms= tot_firms/numtimes ;
totentry = totentry /numtimes ;
totshare0= totshare0/numtimes ;
totshare1= totshare1/numtimes ;
totshare2= totshare2/numtimes ;  % uses 0 for states in which only 1 firm 
totprof1 = totprof1 /numtimes ;
totprof2 = totprof2 /numtimes ;  % uses 0 for states in which only 1 firm 
totprice1= totprice1/numtimes ;
totprice2= totprice2/numtimes ;  % uses 0 for states in which only 1 firm 
totmarkup= totmarkup/numtimes ;
totwstar = totwstar /numtimes ;
totwsame = totwsame /numtimes ;
totLRstay= totLRstay/numtimes ;
totLRx   = totLRx   /numtimes ;
totw2    = totw2    /numtimes ;  % uses 0 for states in which only 1 firm 
totpdiff = totpdiff /numtimes ;
hitkmax  = hitkmax  /numtimes ;

% totexit  = totexit/numtimes;
totv0  = totv0./totexit;
totvx  = totvx./totexit;


disp(['RLGforced = ' num2str(RLGforced) ' times, or ' num2str(RLGforced/numtimes) ' share of simulated periods (last run)']);  % only uses last run since RLGforced,RLGinnov reset each run to zero
disp(['RLGinnov  = ' num2str(RLGinnov)  ' times, or ' num2str(RLGinnov/numtimes)  ' share of simulated periods (last run)']);
disp(' ');

disp(' ');
disp('---------------------------------  Simulated Moments  ---------------------------------');
disp(['Mean disc EU period T : ' mse(discEU_T)]);
disp(['Mean disc CS period T : ' mse(discCS_T)]);
disp(' ');
disp(['Mean consumer surplus : ' mse(consurp)]);
disp(['Mean producer surplus : ' mse(prodsurp)]);
disp(['Mean    total surplus : ' mse(consurp + prodsurp)]);
% disp(' ');
% disp(['Mean CS    (adj. p*)  : ' mse(consurp2)]);         %% ignore these  *2  variables (here and above): from legacy code attempting to approximate equilibrium with unbounded state space
% disp(['Mean PS    (adj. p*)  : ' mse(prodsurp2)]);
% disp(['Mean CS+PS (adj. p*)  : ' mse(consurp2+prodsurp2)]);
disp(' ');
disp(['Mean industry innovate: ' mse(totinnov )]);
disp(['Mean  leader  innovate: ' mse(totinnov1)]);
disp(['Mean  firm 2  innovate: ' mse(totinnov2)]);
disp(['Mean industry invest  : ' mse(totinvest)]);
disp(['Mean  leader  invest  : ' mse(totinves1)]);
disp(['Mean industry invest0 : ' mse(totinves0)]);
disp(['Mean  last t  invest0 : ' mse(totxT0   )]);
disp(['Mean  last t front.inn: ' mse(totinnovT)]);
disp(['Mean freq #firms=  0  : ' mse(tot_firms(:,Params.max_firms+1))]);
for nr=1:Params.max_firms
    disp([ sprintf('Mean freq #firms= %2.0f  : ', nr)  mse(tot_firms(:,nr))]);
end
disp(['Mean leader advantage : ' mse(totw2)]);
disp(['Mean leader advantageT: ' mse(totw2T)]);
disp(['Mean # firms active   : ' mse(totnfirms)]);
disp(['Mean # firms invest0  : ' mse(totinve0n)]);
disp(['Mean # top 2 same w   : ' mse(totwsame )]);
disp(['Mean # firms > wstar  : ' mse(totwstar )]);
disp(['Mean # wstar /nfirms  : ' mse(totwstar./totnfirms)]);
disp(['Mean freq of entry    : ' mse(totentry )]);
disp(['Mean freq lowrung stay: ' mse(totLRstay)]);
disp(['Mean freq lowrung x>0 : ' mse(totLRx   )]);
disp(['Mean leader at kmax   : ' mse(hitkmax)]);
disp(['  max k by any firm   : ' sprintf('%9.0f', maxk)]);
disp(' ');
disp(['Mean share out. good  : ' mse(totshare0)]);
disp(['Mean share leader     : ' mse(totshare1)]);
disp(['Mean share firm 2     : ' mse(totshare2)]);
disp(['Mean share leader @ T : ' mse(totshar1T)]);
disp(['Mean industry markup  : ' mse(totmarkup)]);
disp(['Mean leader   markup  : ' mse(totprice1/mc)]);
disp(['Mean firm 2   markup  : ' mse(totprice2/mc)]);
disp(['Mean leader   profit  : ' mse(totprof1)]);
disp(['Mean firm 2   profit  : ' mse(totprof2)]);
disp(['Mean price(w,.) - p*  : ' mse(totpdiff)]);
disp(' ');
disp(['Mean initial  value   : ' mse(totv0)]);
disp(['Mean realized value   : ' mse(totvx)]);
disp(['Mean initial  v 1st   : ' mse(totv0i)]);
disp(['Mean realized v 1st   : ' mse(totvxi)]);
disp(['  sh 1st firm exits   : ' sprintf(' %12.3f', mean(totlspi<=numtimes))]);
disp('---------------------------------  Done Simulation  ---------------------------------');

%% add whatever you want to return to calling function
simstuff.consurp   = consurp;
simstuff.prodsurp  = prodsurp;
simstuff.totinnov  = totinnov;
simstuff.totinnov1 = totinnov1;
simstuff.totinnov2 = totinnov2;
simstuff.totinnovT = totinnovT;
simstuff.tot_firms = tot_firms;
simstuff.totnfirms = totnfirms;
simstuff.totw2     = totw2;
simstuff.totw2T    = totw2T;
simstuff.totentry  = totentry;
simstuff.hitkmax   = hitkmax;
simstuff.totshare0 = totshare0;
simstuff.totshare1 = totshare1;
simstuff.totshare2 = totshare2;
simstuff.totshar1T = totshar1T;
simstuff.totmarkup = totmarkup;
simstuff.totprof1  = totprof1;
simstuff.totprof2  = totprof2;

end % End of function welf_ma



%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Function encode - takes a weakly descending n-tuple (n =
% Params.max_firms), with min. elt. 0, max. elt. kmax, and encodes it into
% an integer
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function code = encode(binom, Params, ntuple)
code = 1;
for i = 1:Params.max_firms
    digit = ntuple(i);
    code = code+binom(digit+Params.max_firms+1-i,digit+1);
end
end % End of function encode



%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Function acronym - get textual description of types
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function [s,  etype,  itype] = acronym(Par, et, it)
if et == Par.uu_competition
    s = '';
    etype = 'markov perfect nash equilibrium';
elseif et == Par.uu_monopoly
    s = 'm';
    etype = 'monopoly';
elseif et == Par.uu_planner
    s = 's';
    etype = 'social planner';
end

if it == Par.uu_quality
    s = [s 'b'];
    itype = 'quality (differentiated products)';
elseif it == Par.uu_cost
    s = [s 'c'];
    itype = 'marginal cost (homogenous products)';
elseif it == Par.uu_capacity
    s = [s 'p'];
    itype = 'capacity (homogenous products)';
end

if Par.batchjob > 0 
    s = ['bj' num2str(Par.batchjob) '_' s]; 
end
end % End of function acronym



%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Function getint - Get integer between nlow and nhigh with default value
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function n = getint(Par, prompt,  nlow,  nhigh,  default)
n = Par.nocheck-1;
while 1
    n = getnum(prompt, default);
    if (n ~= round(n))
        disp('please re-enter, this parameter must be integer-valued');
        continue;
    end
    if((n<nlow && nlow ~= Par.nocheck) || (n>nhigh && nhigh ~= Par.nocheck))
        disp('please re-enter. Valid range for this parameter is ');
        if(nhigh == Par.nocheck)
            disp(['>' int2str(nlow)]);
        elseif(nlow == Par.nocheck)
            disp(['<' int2str(nhigh)]);
        else
            disp([int2str(nlow) '-' int2str(nhigh)]);
        end
    else
        break;
    end
end
end % End of function getint



%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Helper functions
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% get a number from the user
function ret = getnum(prompt,  default)
ret = input(prompt);
if (length(ret) == 0)
    ret = default;
end
end % End of function getnum

% Implement gauss selif 
function y = selif(x, e)
y = x(e == 1, :);
end

% Printing of statistics on darray
function ret1 = mse(darray)
  ret1 = sprintf(' %12.3f   std: %12.5f   stderr: %12.6f   range: %12.3f  %12.3f ', mean(darray), std(darray), std(darray)/sqrt(size(darray,1)), min(darray), max(darray) );
end
