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
- Use Bellman iteration
- Update the transition rule for N
% 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