Contents

function [nPrime,nPrimeIndex,vFuncs,minN,maxN,Pi,omega,nCstates] = bellman(piF,kappa,beta,phi,markovARG,varargin)

bellman.m

Calculates LIFO equilibrium continuation values using Bellman iteration and plots the results.

Input: piF ~ Producer surplus per customer kappa ~ per-period fixed cost beta ~ annual interest rate phi ~ sunk cost of entry

Output:

Create Markov chain if required.

The results are stored in Pi and omega.

if isfield(markovARG,'Pi')
    Pi = markovARG.Pi;
    omega = markovARG.omega;
    nCstates = markovARG.nCstates;
else

Rip out markovARG parameters if generating

    approximateARG = markovARG.approximateARG;
    omegaCenter = markovARG.omegaCenter;
    omegaWidth = markovARG.omegaWidth;
    omegaStep = markovARG.omegaStep;
    rho = markovARG.rho;
    sigma = markovARG.sigma;
    [Pi,omega,nCstates] = markov(omegaCenter,omegaWidth,omegaStep,rho,sigma,approximateARG);
end
Input argument "markovARG" is undefined.

Error in ==> bellman at 21
if isfield(markovARG,'Pi')

Calculate upper bound on firm numbers.

nCheck=find((max(exp(omega))*piF(1:1:100)./(1:1:100)>=kappa),1,'last');

Solve individual entrants' problems

% Store the individual firms' value functions in a three-dimensional array.
% The final dimension gives the rank of the firm.
% The second dimension gives the number of firms presently active.
% The first dimension indexes the demand states.

vFuncs=zeros(nCstates,nCheck,nCheck);


for R=nCheck:-1:1

    %Initialize the value function.
    v=zeros(nCstates,nCheck-R+1);

    %Initialize the transition rule for N if necessary.
    if R==nCheck;
        nPrime=nCheck*ones(nCstates,1);
    end

    %Create nPrimeIndex, which determines which column of v(:,:) contains the relevant continuation value.
    if R==nCheck
        nPrimeIndex=ones(nCstates,nCheck-R+1);
    else
        nPrimeIndex=zeros(nCstates,nCheck-R+1);
        for i=R:1:nCheck;
            nPrimeIndex=nPrimeIndex+(nPrime==i)*(1+(i-R));
        end
    end

    % Bellman Iteration
    tol=1.0e-7;
    vDistance=1;

    while vDistance>=tol

        %Use nprime to assign the appropriate continuation value.
        Piv=Pi*v;
        Ev=zeros(nCstates,nCheck-R+1);
        for i=1:nCstates;
            Ev(i,:)=Piv(i,nPrimeIndex(i,:));
        end

        vPrime=beta*(Pi*exp(omega')*ones(1,nCheck-R+1)).*piF(nPrime)./nPrime-beta*kappa+beta*Ev;

        %Take the maximum of the expected continuation value and zero.
        vPrime=max(vPrime,0);

        %Measure the distance between the initial and new value functions.
        vDistance=max(max(abs(vPrime-v)));

        %Update the value function.
        v=vPrime;

    end

    %Store the value function
    if R>1; %In this case, we need to pad the matrix of values with NaNs
        vv=[NaN(nCstates,R-1) v];
    else
        vv=v;
    end

    vFuncs(:,:,R)=vv;

    %Update nprime for the problem of a firm with rank R-1

    %Take away one firm if exit is optimal.
    nPrime=nPrime-(v==0);

    %Add a column to nPrime for the initial state N=R-1 (The Matlab editor complains about the speed consequences of
    %this assigment. For a big problem where the speed drag from copying nPrime to a new location in memory is too much
    %to pay, preallocate nPrime to a nCstates x nCheck+1 array and only use the relevant column in the Bellman equation
    %above.)
    nPrime=[(R-1)*ones(nCstates,1) nPrime];

    %Add firms to the new column in the states where entry is optimal.
    for i=R:1:nCheck
        nPrime(:,1)=nPrime(:,1)+(squeeze(vFuncs(:,i,i))>phi(i));
    end

end

Create final version of nPrimeIndex.

nPrimeIndex=zeros(nCstates,nCheck+1);
for i=1:1:nCheck;
    nPrimeIndex=nPrimeIndex+(nPrime==i)*(1+i);
end

Determine the ergodic set for N.

minN=min(min(nPrime));
maxN=max((max(nPrime).*(max(nPrime)>(0:1:nCheck)))); %Highest value of N transited to from a lower value.

%Trim the matrix nPrime to eliminate states outside of the ergodic set.
nPrime=nPrime(:,minN+1:maxN+1);
nPrimeIndex=nPrimeIndex(:,minN+1:maxN+1)-minN+1;

%Trim the value functions to eliminate values
%of N above the ergodic set.
vFuncs=vFuncs(:,1:maxN,1:maxN);

Plot the value functions.

This section creates a raw inelegant plot which can be used to verify that the value functions satisfy the restrictions of theory.

if nargin > 5

    figure(1)
    set(gcf,'Units','inches','Position',[0 0 4 6]);
    for i=1:maxN;
        subplot(maxN,1,i);
        vi=squeeze(vFuncs(:,:,i));
        plot(omega'*ones(1,maxN-i+1),vi(:,i:maxN)); %This line plots the value function's many ``branches'' all at one go.
        hold
        plot(omega',phi(i)*ones(size(omega)));      %This line adds a horizontal line to mark the entry cost.
    end
end

End of File

end