% Originally written in Gauss by Gautam Gowrisankaran in 1993
% modified by Ron Goettler in 2012
% converted to Matlab in 2013
function profit_main(batchnum)

disp('Entering profit_main.m ...');
  
% Build Par struct with general program data
Par.batchjob        = batchnum;
Par.configfile      = ['config_' num2str(Par.batchjob) '.mat'];
Par.small           = 1e-8;
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]);


% Set rlg params and do data tests
Par.Rlg.w0 = 1;
disp(' ');
if Par.Rlg.w0 == 1
    disp('hey: in profit_main.m Rlg.w0 = 1-->lowest w is a place-holder for an open spot, not a firm offering lowest quality');
else
    disp('hey: in profit_main.m Rlg.w0 = 0-->lowest w is actually a firm, not a place-holder for an open spot');
end

if (Params.rlg_y>0 && Params.wstar>0)
    disp(' '); 
    disp ('hey: should not use both rlg_y && wstar'); 
    disp(' '); 
end
if (Params.rlg_y>0 && Params.rlg_y<Params.mc)
    disp(' '); 
    disp ('hey: must have rlg_y>mc else profits = zero'); 
    disp(' '); 
end

disp(' '); 
if (Params.rlg_outgood) 
    disp('hey: outside good is turned on');
else
    disp('hey: outside good is turned off. monopolist price is rlg_y(if rlg_y>0) || uu_rlg_maxp*leader price in max differentiated duopoly');
end
disp(' '); 

if (Params.rlg_sh_cap<.5)
  disp(' '); 
  disp('hey: share cap must be between .5 & 1');  
  disp('Aborting!');
  disp(' ');
  return;
end;

Par.Rlg.wstar   = Params.wstar;
Par.Rlg.maxp    = 9999;
Par.Rlg.fix_sh1 = 0;
Par.Rlg.w       = 0;


% Set kmax based on relevant competition case
kmax = iif(Params.eql_type~=Par.uu_competition,Params.kmax+1,Params.kmax);
kmax = iif(Params.ind_type==Par.uu_capacity && Params.eql_type==Par.uu_planner,kmax+1,kmax);
kmax = iif(Par.Rlg.wstar==0 && Params.eql_type==Par.uu_competition,kmax+1,kmax);


% Set up binomial coefficients for decoding/encoding of n-tuples
Par.Rlg.binom = eye(Params.max_firms+kmax+2);
Par.Rlg.binom = [zeros(Params.max_firms+kmax+2,1) Par.Rlg.binom];
for i=2:(Params.max_firms+kmax+2)
    Par.Rlg.binom(i,2:i) = Par.Rlg.binom(i-1,2:i)+Par.Rlg.binom(i-1,1:i-1);
end


% Start main computation loop
disp('*** computing profit function ***')

for nfirms = 1:Params.max_firms
    % Number of descending n-tuples
    descn = Par.Rlg.binom(nfirms+kmax+1,kmax+2);
    disp(['Number of firms = ' int2str(nfirms) ' ; Total industry structures to compute: ' int2str(descn) ])

    % Call the correct profit function based on current case
    if Params.ind_type == Par.uu_quality
        if (nfirms == 1 && Params.rlg_y == 0)
            Par.Rlg.w = ([kmax; 1])*Params.rlg_wscale+Params.rlg_wshift;
            p = [Params.mc+.6; Params.mc+.5];
            [p,~] = newton(Par,Params,p,'cfunk');  % No need for sigma1 here
            Par.Rlg.maxp = Params.rlg_maxp*p(1);
            disp(['Max differentiated duopoly prices: ',num2str(p'),' ; Max monopoly price is: ',num2str(Params.rlg_maxp),' ; *p_leader = ',num2str(Par.Rlg.maxp)])
        end
        if Params.eql_type == Par.uu_competition
            [profit,agprof,csurplus,share,pmcmarg,concent,price] = cqprofit(Par,Params,nfirms,descn);
        elseif Params.eql_type == Par.uu_monopoly
            [profit,agprof,csurplus,share,pmcmarg,concent,price] = mqprofit(Par,Params,nfirms,descn);
        elseif Params.eql_type == Par.uu_planner
            [profit,agprof,csurplus,share,pmcmarg,concent,price] = sqprofit(Par,Params,nfirms,descn);
        end
    elseif Params.ind_type == Par.uu_cost
        if Params.eql_type == Par.uu_competition
            [profit,agprof,csurplus,share,pmcmarg,concent,price] = ccprofit(Par,Params,nfirms,descn);
        elseif Params.eql_type == Par.uu_monopoly
            [profit,agprof,csurplus,share,pmcmarg,concent,price] = mcprofit(Par,Params,nfirms,descn);
        elseif Params.eql_type == Par.uu_planner
            [profit,agprof,csurplus,share,pmcmarg,concent,price] = scprofit(Par,Params,nfirms,descn);
        end
    elseif Params.ind_type == Par.uu_capacity
        if Params.eql_type == Par.uu_competition
            [profit,agprof,csurplus,share,pmcmarg,concent,price] = cpprofit(Par,Params,nfirms,descn);
        elseif Params.eql_type == Par.uu_monopoly
            [profit,agprof,csurplus,share,pmcmarg,concent,price] = mpprofit(Par,Params,nfirms,descn);
        elseif Params.eql_type == Par.uu_planner
            [profit,agprof,csurplus,share,pmcmarg,concent,price] = spprofit(Par,Params,nfirms,descn);
        end
    end

	% agprof   -> Aggregate profit per industry structure
	% csurplus -> Consumer surplus
	% share    -> Market share of each firm
	% price    -> Price of each firm
	% pmcmarg  -> Price/mc margin, average by sales
	% concent  -> One-firm concentration ratio

    if (nfirms == Params.max_firms  || Par.batchjob == 0)
        file1 = [prefix,'pr.',int2str(nfirms),'f' '.mat'];
        file2 = [prefix,'cons.',int2str(nfirms),'f' '.mat'];
        disp(['generating output-->',file1,',',file2])
        save(file1,'profit');
        save(file2,'agprof','csurplus','share','pmcmarg','concent','price');
    end

    if Par.Rlg.w0 == 1
        Par.Rlg.profit      = profit;
        Par.Rlg.agprof      = agprof;
        Par.Rlg.csurplus    = csurplus;
        Par.Rlg.share       = share;
        Par.Rlg.price       = price;
        Par.Rlg.pmcmarg     = pmcmarg;
        Par.Rlg.concent     = concent;
    end
end

Params.profit_done = 1;
save(Par.configfile, 'Params');
end % End of function profit_main



%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Function cqprofit - profit for quality competition
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function [profit,agprof,csurplus,share,pmcmarg,concent,price] = cqprofit(Par,Params,nfirms,descn)

nagents = iif(Params.eql_type == Par.uu_competition, nfirms, 1);
profit = zeros(descn,nagents);
agprof = zeros(descn,nfirms);
csurplus = zeros(descn,1);
share = zeros(descn,nfirms);
price = zeros(descn,nfirms);
pmcmarg = zeros(descn,1);
concent = zeros(descn,1);

p = (Params.mc+.5)*ones(nfirms,1);

for i=1:descn
    if (mod(i,1000) == 0)
        disp(['computed: ',int2str(i)]);   
    end
    Par.Rlg.w = qdecode(Par,Params,i,nfirms+1);
    if (Par.Rlg.w0 == 1 && Par.Rlg.w(size(Par.Rlg.w,1)) == Params.rlg_wshift)
        price(i,nfirms) = 0;
        if nfirms>1
            k = encode(Par,Params,round((Par.Rlg.w(1:nfirms-1)-Params.rlg_wshift)/Params.rlg_wscale),nfirms);
            profit(i,1:nfirms-1)    = Par.Rlg.profit(k,:);
            agprof(i,1:nfirms-1)    = Par.Rlg.agprof(k,:);
            share(i,1:nfirms-1)     = Par.Rlg.share(k,:);
            price(i,1:nfirms-1)     = Par.Rlg.price(k,:);
            pmcmarg(i)              = Par.Rlg.pmcmarg(k);
            csurplus(i)             = Par.Rlg.csurplus(k);
            concent(i)              = Par.Rlg.concent(k);
        end
    else
        if (Par.Rlg.w(nfirms) == Params.rlg_wshift)
            if (Params.rlg_y>0)
                p = min([(Params.mc+1+.2*Par.Rlg.w'); (Params.rlg_y-.001*Params.rlg_alpha)*ones(1,nfirms)]);
            else
                p(nfirms) = max([Params.mc+.2; Params.mc+1+.2*Par.Rlg.w(nfirms)/Params.rlg_alpha]); 
                for k = nfirms-1:-1:1
                    p(k) = p(k+1)+(Par.Rlg.w(k)-Par.Rlg.w(k+1))/Params.rlg_alpha; 
                end
            end
        end

        if (nfirms == 1)
            if (Params.rlg_y>0)
                if (Params.rlg_outgood == 0)
                    p = Params.rlg_y-Par.small;
                else
                    [p,sigma1] = newton(Par,Params,Params.mc+.5,'cfunk');
                    if (sigma1(1)>Params.rlg_sh_cap)
                        p = Params.rlg_y*(1-exp(log(Params.rlg_sh_cap/((1-Params.rlg_sh_cap)*eg(Par,Params,Par.Rlg.w(1))))/Params.rlg_alpha));
                    end
                end
            else
                if (Params.rlg_outgood == 0)
                    p = Par.Rlg.maxp;
                else
                    [p,sigma1] = newton(Par,Params,Params.mc+.5,'cfunk');
                    if (sigma1(1)>Params.rlg_sh_cap)
                        p = log(Params.rlg_sh_cap/((1-Params.rlg_sh_cap)*eg(Par,Params,Par.Rlg.w(1))))/(-Params.rlg_alpha);
                    end
                end
            end
        else
            [p,sigma1] = newton(Par,Params,p,'cfunk');
            if (sigma1(1)>Params.rlg_sh_cap)
                Par.Rlg.fix_sh1 = 1;
                [p,sigma1] = newton(Par,Params,p(2:size(p,1)),'cfunk');
                p = [0; p];
                Par.Rlg.fix_sh1 = 0;
            end
        end
        
        for j = 1:nfirms-1
            k = (Par.Rlg.w == Par.Rlg.w(j));
            if (sum(k)>1)
                k = find(k == 1)';
                p(k) = mean(p(k))*ones(size(k,1),1);
            end
        end
        [k,sigma1] = cfunk(Par, Params, p);
        profstar = Params.mkt_size*p.*sigma1-Params.mkt_size*Params.mc*sigma1-Params.fixed_cost.*(Par.Rlg.w ~= Params.rlg_wshift);
        profit(i,:) = profstar';
        agprof(i,:) = profstar';
        egw = eg(Par,Params,Par.Rlg.w);
	%   ###########################################   Consumer Surplus issues   #################################################     %
	% RLG_out is just an approximation -- could equivalently have RLG_wshift really high, but might get exp() overflow issues         %
	% TRUE absence of outside good --> utility outside good = -infinity -->  ...                                                      %
	% --> ANY period with monopolist yields -infinity for that period and hence discounted EU is -infinity if EVER get monopolist.    %
	% CS is infinity, or arbitrarily high, if no outside good or if approximating no outside good w/ arbitrarily high RLG_wshift.     %
	%                                                                                                                                 %
	% One "fix" --> redefine utility w/out the industry as 0 if RLG_y = 0 or RLG_alpha * log(RLG_y) if RLG_y > 0.  But CS<0 possible  %
	% Another   --> redefine utility w/out the industry as utility when 1 firm at lowest w.                                           %
	%                                                                                                                                 %
	% OR, in welf_ma.g focus on growth in CS, perhaps as  mean CS_finalperiod / CS_initperiod.  Good for delta=0, but delta>0 ?       %
	% This ratio works well with linear, but less so with log(y-p) since really want to measure Compensating Y needed to get disc.EU  %
	% With log(y-p), ratio of CS_fullrun / (CS_initperiod/(1-beta)) might be appropriate.                                             %
	%                                                                                                                                 %
	% With log(y-p), period CS quickly explodes as  huge y* needed for log(y*) to match  +RLG_wscale * RLG_forced (ie. grid shifts)   %
	%                                                                                                                                 %
	% csurplus below is  M * log(sum(exp( utility ))) /RLG_alpha  including the outside good, which therefore depends on RLG_wshift.  %
	% If RLG_out = 1 utility for out.good, excluding log(y) if present, hardcoded as RLG_wshift-30 to avoid -infinity if monopolist.  %
	%                                                                                                                                 %
	% If RLG_y > 0 then in welf_ma.g will compute the period CS as  M * (exp( csurplus * RLG_alpha )/ RLG_alpha - RLG_y), since       %
	% y* solving   alpha*log( y* ) =  log( sum( exp( utility )))   is income per person needed this period to match EU with industry. %
	%                                                                                                                                 %
	% If RLG_y > 0 cannot compute discounted CS by simply dividing by RLG_alpha, since marginal utility of $ not constant.  Options:  %
	%   1. discounted sum of period CS computed as above:    M * ( exp( csurplus * RLG_alpha )/ RLG_alpha - RLG_y )                   %
	%   2. M*(y*-y)/(1-beta) where  alpha*log(y*)/(1-beta) = disc.sum  beta^t * csurplus * RLG_alpha                                  %
	%                            --> y* = exp( (1-beta)/alpha * sum( beta^t * (csurplus + RLG_forced*RLG_wscale) * RLG_alpha / M )    %
	%   1 is closer to standard/linear case of discounting CS in each period, but 2 is the true compensating variation.               %
        if (Params.rlg_y)
            csurplus(i) = Params.mkt_size*log(exp(log(Params.rlg_y))+sum(exp((1-Params.rlg_outgood)*max([0; 10-Params.rlg_wshift])+log(egw)+Params.rlg_alpha*log(Params.rlg_y-p))))/Params.rlg_alpha;
        else
            csurplus(i) = Params.mkt_size*log(1+sum(exp((1-Params.rlg_outgood)*max([0; 10-Params.rlg_wshift])+log(egw)-Params.rlg_alpha*p)))/Params.rlg_alpha;
        end
        share(i,:) = sigma1';
        price(i,:) = p';
        pmcmarg(i) = sum(p.*sigma1)/Params.mc/sum(sigma1);
        concent(i) = max(sigma1)/sum(sigma1);
    end

    if nfirms<3
        disp(['w: ' sprintf('%8.4f', Par.Rlg.w),'  prices: ' sprintf('%8.4f',price(i,:)) '  mktshares: ' sprintf('%8.5f',share(i,:)) '  sumshares: ' sprintf('%8.5f', sum(share(i,:)')) '  profit: ' sprintf('%8.4f',profit(i,:)) '  cs: ' sprintf('%8.4f', csurplus(i)) ]);
    end
end
end % End of function cqprofit



%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Function mqprofit - profit for quality monopolist
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function [profit,agprof,csurplus,share,pmcmarg,concent,price] = mqprofit(Par,Params,nfirms,descn)

nagents = iif(Params.eql_type == Par.uu_competition, nfirms, 1);
profit = zeros(descn,nagents);
agprof = zeros(descn,nfirms);
csurplus = zeros(descn,1);
share = zeros(descn,nfirms);
price = zeros(descn,nfirms);
pmcmarg = zeros(descn,1);
concent = zeros(descn,1);

p = (Params.mc+.5)*ones(nfirms,1);

if (Params.rlg_alpha ~= 1)
    disp('hey: need to add rlg_alpha to monopoly(cartel) case in profit_main.m... aborting!');
    return;
end

for i=1:descn
    if (mod(i,1000) == 0)
        disp(['computed: ',int2str(i)]);   
    end
    Par.Rlg.w = qdecode(Par,Params,i,nfirms+1);
    Par.Rlg.w = max(([(Par.Rlg.w-Params.rlg_wscale) (Params.rlg_wshift*ones(nfirms,1))])');
	% As zero here represents being out, move everything down by one, 
    % except for zero, as there is no negative efficiency level
    if (nfirms>1 && Par.Rlg.w(nfirms) == Params.rlg_wshift)
        if (Params.rlg_y>0) 
            p = min(([Params.mc+.5*(Par.Rlg.w+1) (Params.rlg_y-.001)*ones(nfirms,1)])');
        else
            p = (Params.mc+1+.2*Par.Rlg.w);
        end
    end
    [p,sigma1] = newton(Par,Params,p,'mfunk');
    profstar = Params.mkt_size*p.*sigma1-Params.mkt_size*Params.mc*sigma1;
    profit(i) = sum(profstar-Params.fixed_cost.*(Par.Rlg.w ~= Params.rlg_wshift));
    agprof(i,:) = profstar';
    disp('hey: need to check whether fixed costs are correctly being deducted only for active plants in mqprofit()in profit_main.m')

    egw = eg(Par,Param,Par.Rlg.w);
    if Params.rlg_y>0
        csurplus(i) = Params.mkt_size*log(Params.rlg_outgood+sum(exp(Par.Rlg.w+log(Params.rlg_y-p)-log(Params.rlg_y))));
    else
        csurplus(i) = Params.mkt_size*log(Params.rlg_outgood+sum(exp(log(egw)-p)));
    end
    share(i,:) = sigma1';
    price(i,:) = p';
    pmcmarg(i) = sum(p.*sigma1)/Params.mc/sum(sigma1);
    concent(i) = max(sigma1)/sum(sigma1);
    if nfirms<3
        disp(['w: ',num2str(Par.Rlg.w'),' ; prices: ',num2str(price(i,:)),' ; mktshares: ',num2str(share(i,:)),' ; sumshares: ',num2str(sum(share(i,:)')),' ; profit: ',num2str(profit(i,:)),' ; cs: ',num2str(csurplus(i))]); 
    end
end
end % End of function mqprofit



%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Function sqprofit - profit for quality planner
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function [profit,agprof,csurplus,share,pmcmarg,concent,price] = sqprofit(Par,Params,nfirms,descn)

nagents = iif(Params.eql_type == Par.uu_competition, nfirms, 1);
profit = zeros(descn,nagents);
agprof = zeros(descn,nfirms);
csurplus = zeros(descn,1);
share = zeros(descn,nfirms);
price = zeros(descn,nfirms);
pmcmarg = zeros(descn,1);
concent = zeros(descn,1);

p = (Params.mc+.5)*ones(nfirms,1);

for i=1:descn
    if (mod(i,1000) == 0)
        disp(['computed: ',int2str(i)]);   
    end
    Par.Rlg.w = qdecode(Par,Params,i,nfirms+1);
    Par.Rlg.w = max(([(Par.Rlg.w-Params.rlg_wscale) (Params.rlg_wshift*ones(nfirms,1))])');
	% As zero here represents being out, move everything down by one, 
    % except for zero, as there is no negative efficiency level
    p = Params.mc*ones(nfirms,1);
    egw = eg(Par,Params,Par.Rlg.w);
    if Params.rlg_y>0 
        egwp = exp(Par.Rlg.w+log(Params.rlg_y-p)-log(Params.rlg_y));
    else
        egwp = exp(log(egw)-p);
    end
    n = egwp;
    sigma1 = n./(1+sum(n));
    % Consumer surplus - plant fixed costs
    profit(i) = Params.mkt_size*log(Params.rlg_outgood+sum(egwp))-Params.fixed_cost*sum(Par.Rlg.w ~= Params.rlg_wshift);
    disp('hey:  need to check whether fixed costs are correctly being deducted only for active plants in sqprofit()in profit_main.m')
    share(i,:) = sigma1';
    price(i,:) = p';
    pmcmarg(i) = sum(p.*sigma1)/Params.mc/sum(sigma1);
    concent(i) = max(sigma1)/sum(sigma1);
    if nfirms<3
        disp(['w: ',num2str(Par.Rlg.w'),' ; prices: ',num2str(price(i,:)),' ; mktshares: ',num2str(share(i,:)),' ; sumshares: ',num2str(sum(share(i,:)')),' ; profit: ',num2str(profit(i,:)),' ; cs: ',num2str(csurplus(i))]); 
    end
end
end % End of function sqprofit



%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% The following three functions compute static profit function for 
% "Cournot" model (homogenous products, investment in marginal cost) for 
% different types of industry conduct (competition, multi-plant monopolist,
% multi-plant social planner).
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Function ccprofit - profit for cost competition
% Solve for equilibrium with n firms; reduce n until all firms want to 
% produce quantity > 0
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function [profit,agprof,csurplus,share,pmcmarg,concent,price] = ccprofit(Par,Params,nfirms,descn)

nagents = iif(Params.eql_type == Par.uu_competition, nfirms, 1);
profit = zeros(descn,nagents);
agprof = zeros(descn,nfirms);
csurplus = zeros(descn,1);
share = zeros(descn,nfirms);
price = zeros(descn,nfirms);
pmcmarg = zeros(descn,1);
concent = zeros(descn,1);

for i=1:descn
    if (mod(i,1000) == 0)
        disp(['computed: ',int2str(i)]);   
    end
    w = cdecode(Par,Params,i,nfirms+1);
    theta = Params.gamma*exp(-w(1)); % Marginal cost
    n = nfirms;
    p = (Params.intercept+sum(theta(1:n)))/(n+1);
    while ~(p-theta(n) >= 0 || n == 1)
        n = n-1;
        p = (Params.intercept+sum(theta(1:n)))/(n+1);
    end;
    q = zeros(nfirms,1);
    if (p-theta(n)>0)
        q(1:n) = p-theta(1:n);
    end
    quan = q;
    pstar = Params.intercept-sum(quan); %  Equilibrium price
    profstar = (pstar>theta).*(pstar-theta).*quan-Params.fixed_cost; % Equilibrium profits
    profit(i,:) = profstar';
    csurplus(i) = 0.5*sum(quan)*sum(quan);
    agprof(i,:) = profstar';
    share(i,:) = quan';
    if sum(quan)>0
        pmcmarg(i) = pstar/(sum(theta.*quan))*sum(quan);
        concent(i) = max(quan)/sum(quan);
    else
        pmcmarg(i) = 1;
    end
end
end % End of function ccprofit



%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Function mcprofit - profit for cost monopoly
% The monopolist will always choose to produce everything from the lowest 
% priced firm, so it acts like a 1-plant firm for the static profits
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function [profit,agprof,csurplus,share,pmcmarg,concent,price] = mcprofit(Par,Params,nfirms,descn)

nagents = iif(Params.eql_type == Par.uu_competition, nfirms, 1);
profit = zeros(descn,nagents);
agprof = zeros(descn,nfirms);
csurplus = zeros(descn,1);
share = zeros(descn,nfirms);
price = zeros(descn,nfirms);
pmcmarg = zeros(descn,1);
concent = zeros(descn,1);

for i=1:descn
    if (mod(i,1000) == 0)
        disp(['computed: ',int2str(i)]);   
    end
    w = cdecode(Par,Params,i,nfirms+1);
    agprof(i,:) = -Params.fixed_cost*((w>-4)');
    numin = sum(w>-4);
    w = max(([(w-1) (-4*ones(nfirms,1))])');
	% As zero here represents being out, move everything down by one, 
    % except for zero, as there is no negative efficiency level
    theta = Params.gamma*exp(-w(1)); % Marginal cost
    pstar = 0.5*(Params.intercept+theta); % One-plant monopolist price
    quan = iif(pstar>theta,pstar-theta,0);
    profstar = quan*(pstar-theta)-Params.fixed_cost*numin; % Monopolist profits
    profit(i) = profstar;
    agprof(i,1) = agprof(i,1)+profstar;
    csurplus(i) = 0.5*quan*quan;
    pmcmarg(i) = pstar/theta;
    concent(i) = 1;
    if nfirms>1
        quan = [quan; zeros(nfirms-1,1)];
    end
    share(i,:) = quan';
end
end % End of function mcprofit



%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Function scprofit - profit for cost planner
% The social planner will always choose to produce everything from the 
% lowest priced firm, so it acts like a 1-plant firm for the static profits
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function [profit,agprof,csurplus,share,pmcmarg,concent,price] = scprofit(Par,Params,nfirms,descn)

nagents = iif(Params.eql_type == Par.uu_competition, nfirms, 1);
profit = zeros(descn,nagents);
agprof = zeros(descn,nfirms);
csurplus = zeros(descn,1);
share = zeros(descn,nfirms);
price = zeros(descn,nfirms);
pmcmarg = zeros(descn,1);
concent = zeros(descn,1);

for i=1:descn
    if (mod(i,1000) == 0)
        disp(['computed: ',int2str(i)]);   
    end
    w = cdecode(Par,Params,i,nfirms+1);
    agprof(i,:) = -Params.fixed_cost*((w>-4)');
    numin = sum(w>-4);
    w = max(([(w-1) (-4*ones(nfirms,1))])');
	% As zero here represents being out, move everything down by one, 
    % except for zero, as there is no negative efficiency level
    theta = Params.gamma*exp(-w(1)); % Marginal cost
    pstar = theta; % Set price = mc, for social planner solution
    quan = iif(Params.intercept>theta,Params.intercept-pstar,0);
    profstar = 0.5*(Params.intercept-theta)*quan; % Consumer surplus
    profit(i) = profstar-Params.fixed_cost*numin; % Producer surplus
    csurplus(i) = profstar;
    pmcmarg(i) = 1;
    concent(i) = 1;
    if nfirms>1
        quan = [quan; zeros(nfirms-1,1)];
    end
    share(i,:) = quan';
end
end % End of function scprofit



%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% The following three functions compute static profit function for 
% "Capacity" model (homogenous products, investment in capacity) for 
% different types of industry conduct (competition, multi-plant monopolist,
% multi-plant social planner).
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Function cpprofit - profit for capacity competition
% Solve for equilibrium with n firms; reduce n until all firms want to 
% produce quantity > 0
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function [profit,agprof,csurplus,share,pmcmarg,concent,price] = cpprofit(Par,Params,nfirms,descn)

nagents = iif(Params.eql_type == Par.uu_competition, nfirms, 1);
profit = zeros(descn,nagents);
agprof = zeros(descn,nfirms);
csurplus = zeros(descn,1);
share = zeros(descn,nfirms);
price = zeros(descn,nfirms);
pmcmarg = zeros(descn,1);
concent = zeros(descn,1);

for i=1:descn
    if (mod(i,1000) == 0)
        disp(['computed: ',int2str(i)]);   
    end
    w = pdecode(Par,Params,i,nfirms+1);
    cap = w*Params.tau;
    n = nfirms;
    sub = 0;
    q = Params.intercept/(n+1)*ones(n,1);
    while q(n)>cap(n)
        q(n) = cap(n);
        sub = sub+q(n);
        n = n-1;
        if n == 0 
            break; 
        end
        q(1:n) = (d-sub)/(n+1)*ones(n,1);
    end;
    quan = q;
    pstar = Params.intercept-sum(quan);
    profstar = pstar.*quan;
    profit(i,:) = profstar';
    csurplus(i) = 0.5*sum(quan)*sum(quan);
    agprof(i,:) = profstar';
    share(i,:) = quan';
    if sum(quan)>0
        pmcmarg(i) = (pstar+Params.mc)/Params.mc;
        concent(i) = max(quan)/sum(quan);
    else
        pmcmarg(i) = 1;
    end
end
end % End of function cpprofit



%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Function mpprofit - profit for capacity monopolist
% Monopolist produces quantity D/2, and charges price on demand curve, 
% unless total capacity is less than D/2, in which case quantity = capacity
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function [profit,agprof,csurplus,share,pmcmarg,concent,price] = mpprofit(Par,Params,nfirms,descn)

nagents = iif(Params.eql_type == Par.uu_competition, nfirms, 1);
profit = zeros(descn,nagents);
agprof = zeros(descn,nfirms);
csurplus = zeros(descn,1);
share = zeros(descn,nfirms);
price = zeros(descn,nfirms);
pmcmarg = zeros(descn,1);
concent = zeros(descn,1);

for i=1:descn
    if (mod(i,1000) == 0)
        disp(['computed: ',int2str(i)]);   
    end
    w = pdecode(Par,Params,i,nfirms+1);
    w = max(([(w-1) zeros(nfirms,1)])');
	% As zero here represents being out, move everything down by one, 
    % except for zero, as there is no negative efficiency level
    cap = w*Params.tau;
    % Solve for total production, and then allocate the production among 
    % the firms putting the maximum possible production on the first firms
    q = zeros(nfirms,1);
    totprod = min([sum(cap); 0.5*Params.intercept]);
    j = 1;
    while(totprod>0) && (j <= nfirms)
        q(j) = min([totprod; cap(j)]);
        totprod = totprod-q(j);
        j = j+1;
    end
    quan = q;
    pstar = Params.intercept-sum(quan);
    profit(i) = sum(quan)*pstar;
    agprof(i,:) = pstar*quan';
    pmcmarg(i) = (pstar+Params.mc)/Params.mc;
    csurplus(i) = 0.5*sum(quan)*sum(quan);
    if sum(quan)>0
        concent(i) = quan(1)/sum(quan);
        share(i,:) = quan';
    end
end
end % End of function mpprofit



%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Function spprofit - profit for capacity planner
% Social planner produces quantity D, and charges price on demand curve, 
% unless total capacity is less than D, in which case quantity = capacity 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function [profit,agprof,csurplus,share,pmcmarg,concent,price] = spprofit(Par,Params,nfirms,descn)

nagents = iif(Params.eql_type == Par.uu_competition, nfirms, 1);
profit = zeros(descn,nagents);
agprof = zeros(descn,nfirms);
csurplus = zeros(descn,1);
share = zeros(descn,nfirms);
price = zeros(descn,nfirms);
pmcmarg = zeros(descn,1);
concent = zeros(descn,1);

for i=1:descn
    if (mod(i,1000) == 0)
        disp(['computed: ',int2str(i)]);   
    end
    w = pdecode(Par,Params,i,nfirms+1);
    w = max(([(w-1) zeros(nfirms,1)])');
	% As zero here represents being out, move everything down by one, 
    % except for zero, as there is no negative efficiency level
    cap = w*Params.tau;
    % Solve for total production, and then allocate the production among 
    % the firms putting the maximum possible production on the first firms
    q = zeros(nfirms,1);
    totprod = min([sum(cap); Params.intercept]);
    j = 1;
    while(totprod>0) && (j <= nfirms)
        q(j) = min([totprod; cap(j)]);
        totprod = totprod-q(j);
        j = j+1;
    end
    quan = q;
    pstar = Params.intercept-sum(quan);
    agprof(i,:) = pstar*quan';
    csurplus(i) = 0.5*sum(quan)*sum(quan);
    profit(i) = sum(agprof(i,:)')+csurplus(i);
    pmcmarg(i) = (pstar+Params.mc)/Params.mc;
    if sum(quan)>0
        concent(i) = quan(1)/sum(quan);
        share(i,:) = quan';
    end
end
end % End of function spprofit



%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Decode functions - take a previously encoded number, and decodes it
% into a weakly descending n-tuple (n = nfirms - 1)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function ntuple = qdecode(Par, Params, code, nfirms)
ntuple = decode(Par, Params, code, nfirms);
% Now convert to format of starting at wshift, and jumping by wscale
ntuple = (ntuple.*Params.rlg_wscale+Params.rlg_wshift);
end % End of function qdecode

function ntuple = cdecode(Par, Params, code, nfirms)
ntuple = decode(Par, Params, code, nfirms);
% Now convert to format of starting at -4, and jumping by 1's
ntuple = ntuple-4;
end % End of function cdecode

function ntuple = pdecode(Par, Params, code, nfirms)
ntuple = decode(Par, Params, code, nfirms);
end % End of function pdecode

function ntuple = decode(Par, Params, code, nfirms)
code = code-1;
ntuple = zeros(nfirms-1,1);
for i=1:(nfirms-1)
    digit = 0;
    while Par.Rlg.binom(digit+nfirms-i+1,digit+2) <= code
        digit = digit+1;
    end
    ntuple(i) = digit;
    code = code-Par.Rlg.binom(digit+nfirms-i,digit+1);
end
end % End of function decode

function code = encode(Par, Params, ntuple, nfirms)
code = 1;
for i = 1:(nfirms-1)
    digit = ntuple(i);
    code = code+Par.Rlg.binom(digit+nfirms-i,digit+1);
end
end % End of function encode



%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Function newton - Newton-Raphson search to find the root of the function
% objfunk. LB <= p <= UB is enforced.
% Note:  objfunk is a function name and must be passed as a string
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function [p,sigma1] = newton(Par, Params, p, objfunk)

% Set newton procedure constants
epsilon     = .0001;   
tol         = 1e-6;   
np          = size(p,1);   
minstep     = 1;            % stepsize is random uniform between [minstep,1]
trymax      = 200*np;       % trymax = iterations before restart search at rand()

% Set upper and lower bounds based on competition type
if Params.ind_type == Par.uu_quality
    lb = zeros(np,1)+Params.mc;
    if Params.rlg_y>0
        ub = zeros(np,1)+Params.rlg_y-Par.small;
    else
        ub = zeros(np,1)+Par.Rlg.maxp;  % Par.Rlg.maxp set as Params.rlg_maxp*price of maximal leader in duopoly
    end
else
    disp('hey: in newton() you might need to add ub & lb definitions for cournot versions of static profit game')
    ub = zeros(np,1)+99999;
    lb = zeros(np,1);
end

% Start the main newton loop
maxiter = trymax*20;
iter = 0;
deriv = zeros(np,np);
while(iter<maxiter)
    iter = iter+1;
    [x,sigma1] = feval(objfunk,Par,Params,p);
    
    % FOC achieved
    if max(abs(x))<tol
        break; 
    end
    
    % Calculate derivative matrix
    for i = 1:np
        dp = p;
        if (ub(i)-p(i)<0.0001) 
            dp(i) = p(i)-epsilon;
        else
            dp(i) = p(i)+epsilon;
        end
        [dx,sigma1] = feval(objfunk,Par,Params,dp);
        deriv(:,i) = (dx-x)/(dp(i)-p(i));
    end

    scalerr = 0;
    try pnew = (inv(deriv));
    catch scalerr;
    end

    if ((mod(iter,trymax)<4 || mod(iter,trymax)>(trymax-4)) && iter>10)
        disp(['newton() at iter ' int2str(iter) ' with w = ' num2str(Par.Rlg.w') ' p = ' num2str(p') ' foc = ' num2str(x') ' deriv = ' num2str(deriv(:)') ' sigma1 = ' num2str(sigma1')]);
    end
    
    if(scalerr ~= 0 || mod(iter,trymax) == 0)
        if (scalerr~=0)
            iter = trymax*(1+floor(iter/trymax));
            disp(['restarting newton() at rand() since singular deriv at iter ' int2str(iter) ' with w = ' int2str(Par.Rlg.w') ' p = ' num2str(p') ' foc = ' num2str(x') ' deriv = ' num2str(deriv(:)') ' sigma1 = ' num2str(sigma1') ]);
        else
            disp(['restarting newton() at rand() since not converging at iter ' int2str(iter) ' with w = ' int2str(Par.Rlg.w') ' p = ' num2str(p') ' foc = ' num2str(x') ]);
        end
        pnew = lb+sortrows(rand(np,1),-1);
        minstep = max([.1; minstep-.1])';
    else
        if mod(iter,trymax)<5
            pnew = pnew*((.1+(mod(iter,trymax))/10)*rand(1,1))*x;
        else
            pnew = pnew*(minstep+(1-minstep)*rand(1,1))*x;
        end
        pnew = p-pnew.*(10/max([10; abs(pnew)])');  % never change price by more than 10 in one step
        if (size(p,1)>1 && sigma1(1)<.001 && p(1)>Params.mc+10)
            pnew(1) = max(pnew(1)-[1; pnew(2)]+.1)';
        end
        if (size(p,1)>1 && mod(iter,trymax)<5 && Params.rlg_y>0)
            pnew(1) = min([pnew(1); (Params.rlg_y-.5+(mod(iter,trymax))/10)])';
        end
    end
    
    pnew = min([ub'; pnew'])';
    pnew = max([lb'; pnew'])';
    p = pnew;
end

if max(abs(x))>tol
    disp(['newton() failed: foc = ' num2str(x') ' ; Returning p = ' num2str(p')])
end
end % End of function newton



%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Function mfunk - used for quality-investment, monopoly model
% Calculate the profit derivative wrt price for the multi-plant monopolist.
% The derivative is as in (18) in PGM paper. However, we can simplify this
% expression as follows:
% first, divide through by sigma(n).
% left with: FOCn = -(pn-mc)(1-sigman) + 1 + sum(k<>n)(sigmak(pk-mc))
% second, define A = sum( sigmak(pk-mc) ).
% left with: FOCn = -(pn-mc)(1-sigman) + 1 + A - sigman(pn-mc).
%   ==> FOCn = -(pn-mc) + 1 + A
%   ==> FOC = -(p-mc) + 1 + A.
%
%  when RLG_y>0, the simplification is .
%  FOCn = -(pn-mc)(1-sigman)/(RLG_y-p) + 1 + sum(k<>n)(sigmak(pk-mc))/(RLG_y-p)
%  multiply by RLG_y-p and then use same A as above when RLG_y=0 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function [ret1,sigma1] = mfunk(Par, Params, p)
if (Params.rlg_y>0)
    n = exp(Par.Rlg.w+log(Params.rlg_y-p)-log(Params.rlg_y));
elseif (Par.Rlg.wstar>0)
    n = eg(Par,Params,Par.Rlg.w).*exp(-p);
else
    n = exp(Par.Rlg.w-p);
end
sigma1 = n./(1.0+sum(n));
a = sum(sigma1.*(p-Params.mc));
if (Params.rlg_y>0)
    ret1 = -(p-Params.mc) + (Params.rlg_y-p) + a;
else
    ret1 = -(p-Params.mc) + 1 + a;
end
end % End of function mfunk



%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Function cfunk - used for quality competition profit function, net of
% fixed costs
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function [ret0,sigma1] = cfunk(Par, Params, p)
if Par.Rlg.fix_sh1
    p = [1; p];
end
if (Params.rlg_y>0) 
    n = eg(Par,Params,Par.Rlg.w).*exp(Params.rlg_alpha*(log(Params.rlg_y-p)-log(Params.rlg_y)));
else
    n = eg(Par,Params,Par.Rlg.w).*exp(-Params.rlg_alpha*p);
end

if Par.Rlg.fix_sh1
    if (Params.rlg_y>0) 
        fixp1 = Params.rlg_y*(1-exp(log((Params.rlg_outgood+sum(n(2:size(n,1))))*Params.rlg_sh_cap/((1-Params.rlg_sh_cap)*eg(Par,Params,Par.Rlg.w(1))))/Params.rlg_alpha));
    else
        fixp1 = log((Params.rlg_outgood+sum(n(2:size(n,1))))*Params.rlg_sh_cap/((1-Params.rlg_sh_cap)*eg(Par,Params,Par.Rlg.w(1))))/(-Params.rlg_alpha);
    end
    p(1) = fixp1;
    n(1) = eg(Par,Params,Par.Rlg.w(1)).*exp(-Params.rlg_alpha*p(1));
end
sigma1 = n./(Params.rlg_outgood+sum(n));
if (Params.rlg_y>0)
    foc = -Params.rlg_alpha*(p-Params.mc).*(1-sigma1)./(Params.rlg_y-p)+1;
else
    foc = -Params.rlg_alpha*(p-Params.mc).*(1-sigma1)+1;
end
ret0 = 	foc(Par.Rlg.fix_sh1+1:size(foc,1));
end % End of function cfunk



%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Function eg - Calculates e^g(w), used for quality competition profit
% function
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function   wret = eg(Par, Params, w)
% checking Par.Rlg.wstar is probably redundant now since checking it 
% before calling eg(). To avoid overflows with large w and large p:  
% exp(w).*exp(-p) less accurate than exp(w-p)
wret = exp(w);
if Par.Rlg.wstar>0
    for i = 1:size(w,1)
        if w(i)>Params.wstar
            % w is shifted, scaled by wscale, wshift
            wret(i) = exp(Params.wstar)*(2.0-exp(-(w(i)-Params.wstar)));
        else
            break; % w sorted High to Low, so no need to check other i
        end
    end
end
end % End of function eg



%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% 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



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

% C-style numeric if else
function ret = iif(a, b, c)
if (a)
    ret =    b;
else
    ret =    c;
end
end % End of function iif

