%Solving Hopenhayn's (1992) firm entry exit model in continuous time.%By SAEED SHAKER-AKHTEKHANE.%This is the 'COMBINED' version where the processes for shock and age are%combined in one matrix to solve for only one Kolmogorov Forward equation.%I have tried to be as consistent as possible with the notation used by%Benjamin Moll in the HACT project. (Have used same sparse matrix method)
clear;close all;warning off;tic;%% SET VALUES%------------Basic parameters and functions--------------------------------gamma1 = 2; % Households CRRA utility with parameter gammaalpha = 0.5; % alpha = 0.35; % labor share in productionqfun = @(shock1,L) exp(shock1).*(L.^(alpha)); % Productionwfun = @(pp,shock1,L) 1; %wagepfun = @(Q1) Q1^(1/(1-gamma1));nfun = @(pp,ww,shock1) ((shock1.^-1).*(ww/(pp*alpha))).^(1/(alpha-1));
zmean = .5; % mean O-U process.sig2 = (0.1)^2; % sigma^2 O-UCorr = exp(-0.3); % persistence -log(Corr) O-Urho = 0.05; % discount rate
N=.5; % initial aggregate labor.Q=qfun(zmean,N); % initial aggregate quantity.relax = .5; % relaxation parameterJ=40; % number of z gridpointszmin = 0; zmax = 1; %shock rangeI=100; % number of age gridpointsamin = 0; amax = 100; %age range
maxit = 100; %maximum number of iterations in the HJB loopmaxitNQ = 60; %maximum number of iterations in the N and Q loopcrit = 1e-6; %accuracy: HJB loopcritQ = 1e-6; %accuracy: Q loopcritN = 1e-6; %accuracy: N loopDelta = 10e5; %delta for HJBDeltag = 10e10; %delta for KFE
theta11 = -log(Corr); % theta: O-UVar = sig2/(2*theta11);
a = linspace(amin,amax,I)'; % age vectorda = a(2)-a(1);da2 = da^2;z = linspace(zmin,zmax,J); % shock vectordz = z(2)-z(1);dz2 = dz^2;aa = repmat(a,1,J); %a*ones(1,J);zz = repmat(z,I,1); %ones(I,1)*z;
mu = theta11*(zmean - z); %Drifts2 = sig2.*ones(1,J); %Variance
%--------------------------------------------------------------------------%Fixed cost (paid every period) and entry cost (paid only once) are%important parameters and main drivers of exit/entry
cfAll = 1.2; % All values for fixed cost. e.g [0:.5:15];ceAll = 2.2; % All values for entry cost. e.g [0:.2:5];
% OR could be creating some logspace for ce (Entry Cost)% % % ce1=.01; ce2=1.5;% % % ce0=ce2/2;% % % H1=ce0-logspace(log10(ce0),log10(ce1),10);% % % H2=ce2 - sort(H1,'descend');% % % H2=H2(2:end);% % % ceAll = [H1,H2];%%--------------------------------------------------------------------------
N_cf = length(cfAll);N_ce = length(ceAll);
XR = zeros(N_cf,N_ce); % Exit/Entry rate matrixTN=zeros(N_cf,N_ce); % Agg. LaborTQ=zeros(N_cf,N_ce); % Agg. QuantityTp=zeros(N_cf,N_ce); % PriceTw=zeros(N_cf,N_ce); % Wage
gz0 = ones(I*J,1);
%% loops on different values for fixed and entry cost
for cfi = 1:N_cf % fixed cost loop for cei = 1:N_ce % entry cost loop cf = cfAll(cfi); %Fixed cost of production - each period c_ex = 0; %Exit cost/fixed payment c_en = ceAll(cei); %Entry cost
% creating diagonal elements of matrix Aswitch (for HJB: from z) yy = - s2/dz2 - mu/dz; zeta = s2/(2*dz2); chi = mu/dz + s2/(2*dz2); % creating diagonal elements of matrix Aswitchg (for KFE: from z) yyg = - s2/dz2 + mu/dz + theta11*ones(1,J); zetag = s2/(2*dz2); chig = - mu/dz + s2/(2*dz2); %creating both sparse matrices updiag=zeros(I,1); updiagg=zeros(I,1); for j=1:J updiag=[updiag;repmat(zeta(j),I,1)]; updiagg=[updiagg;repmat(zetag(j),I,1)]; end centdiag=repmat(chi(1)+yy(1),I,1); centdiagg=repmat(chig(1)+yyg(1),I,1); for j=2:J-1 centdiag=[centdiag;repmat(yy(j),I,1)]; centdiagg=[centdiagg;repmat(yyg(j),I,1)]; end centdiag=[centdiag;repmat(yy(J)+zeta(J),I,1)]; centdiagg=[centdiagg;repmat(yyg(J)+zetag(J),I,1)]; lowdiag=repmat(chi(2),I,1); lowdiagg=repmat(chig(2),I,1); for j=3:J lowdiag=[lowdiag;repmat(chi(j),I,1)]; lowdiagg=[lowdiagg;repmat(chig(j),I,1)]; end Aswitch=spdiags(centdiag,0,I*J,I*J)+... spdiags(lowdiag,-I,I*J,I*J)+spdiags(updiag,I,I*J,I*J); Aswitchg=spdiags(centdiagg,0,I*J,I*J)+... spdiags(lowdiagg,-I,I*J,I*J)+spdiags(updiagg,I,I*J,I*J);
% creating diagonal elements of matrix AA (for HJB: from a) X = (1/da)*ones(I,J); Z = zeros(I,J); Y = (-1/da)*ones(I,J); % creating matrix AA updiag=[0]; for j=1:J updiag=[updiag;Z(1:I-1,j);0]; end centdiag=reshape(Y,I*J,1); lowdiag=X(2:I,1); for j=2:J lowdiag=[lowdiag;0;X(2:I,j)]; end AA=spdiags(centdiag,0,I*J,I*J)+... spdiags([updiag;0],1,I*J,I*J)+spdiags([lowdiag;0],-1,I*J,I*J); % Adding parts of HJB from a and z to get matrix A A = AA + Aswitch;
% creating diagonal elements of matrix AA_g (KFE: from a (NO EXIT)) Z_g = (1/(da))*aa; Z_g(1,:)=(1/(da))*ones(1,J); X_g = zeros(I,J); Y_g = (-1/(da))*aa; Y_g(1,:)= -Z_g(1,:);
% creating diagonal elements of matrix AA_g_E (KFE: from a (EXIT)) Z_g_E = zeros(I,J); X_g_E = (1/da2-1/da)*aa; X_g_E(1,:)=(1/(da))*ones(1,J); Y_g_E = (-1/da2)*aa; Y_g_E(1,:)= -X_g_E(1,:);
% creating matrices AA_g and AA_g_E updiag_g=[0]; updiag_g_E=[0]; for j=1:J updiag_g=[updiag_g;Z_g(1:I-1,j);0]; updiag_g_E=[updiag_g_E;Z_g_E(1:I-1,j);0]; end centdiag_g=reshape(Y_g,I*J,1); centdiag_g_E=reshape(Y_g_E,I*J,1); lowdiag_g=X_g(2:I,1); lowdiag_g_E=X_g_E(2:I,1); for j=2:J lowdiag_g=[lowdiag_g;0;X_g(2:I,j)]; lowdiag_g_E=[lowdiag_g_E;0;X_g_E(2:I,j)]; end AA_g=spdiags(centdiag_g,0,I*J,I*J)+... spdiags([updiag_g;0],1,I*J,I*J)+... spdiags([lowdiag_g;0],-1,I*J,I*J); AA_g_E=spdiags(centdiag_g_E,0,I*J,I*J)+... spdiags([updiag_g_E;0],1,I*J,I*J)+... spdiags([lowdiag_g_E;0],-1,I*J,I*J); % Adding the parts of KFE from a and z with NO Exit A_g = AA_g + Aswitchg; % Adding the parts of KFE from a and z with Exit A_g_E = AA_g_E + Aswitchg;
%% Initial Values p = 1; w = wfun(p,zmean,N); n = nfun(p,w,zz); q = qfun(zz,n); pi0 = p*q-w*n-w*cf; %profits v0r = zeros(I,J); %value for ivr = 1:I v0r(ivr,:) = exp(-rho*a(1:ivr)')*pi0(1:ivr,:); end v0 = p*q-w*n-w*cf; V = v0; dist = zeros(1,maxit); % to check value function convergence
flagN=0; % This is to verify convergence of N flagQ=0; % This is to verify convergence of Q % Main loop ------------------------------------------------------- for iterNQ=1:maxitNQ disp(iterNQ); %% solving HJBVI using LCP solver for ii=1:maxit n = nfun(p,w,zz); q = qfun(zz,n); pi1 = p*q-w*n-w*cf; pi_stacked = reshape(pi1,I*J,1); V_stacked = reshape(V,I*J,1); S= - c_ex*ones(I,J) - c_en*exp(rho*aa); S_st = reshape(S,I*J,1); B = rho*speye(I*J) - A; qq = -pi_stacked + B*S_st; zzv0 = V_stacked - S_st; lb = zeros(I*J,1); ub = Inf*ones(I*J,1); zzv = LCP(B,qq,lb,ub,zzv0,0); LCP_error = max(abs(zzv.*(B*zzv + qq))); if LCP_error > crit disp('LCP not solved') end V_LCP = zzv+S_st; error = zzv.*(B*zzv + qq); V_LCP1 = reshape(V_LCP,I,J); Vchange = V_LCP1 - V; V = V_LCP1; dist(ii) = max(max(abs(Vchange))); if dist(ii)<crit break end end % Finding Exit rate ------------------------------------------- [mlocX,adjX] = find(abs(zzv)==0); [mlocN,adjN] = find(abs(zzv)>=crit); mlocX=mlocX'; mlocN=mlocN'; Xprob0=length(mlocX)/length(zzv); Xprob=sum(gz0(mlocX))/sum(gz0); %Xprob = sum(gz0(1:floor(Xprob0*J)))/sum(gz0); %% Solving KFE ------------------------------------------------ A_gX = (1-Xprob)*A_g + (Xprob)*A_g_E; %Combine exit and No exit AT = (1/Deltag)*speye(I*J) - A_gX'; g=10*ones(I,J); % distribution over z and a gg=reshape(g,I*J,1); for ii = 1:maxit b = (1/Deltag)*gg; gg_prime = AT\b; gg_error = max(abs(gg_prime - gg)); gg = gg_prime; if gg_error < crit; break; end end g_sum = sum(gg); if g_sum ~= 0; gg = gg./g_sum; end g = reshape(gg,I,J); gz0=gg; %% Update N and Q n_st=reshape(n,I*J,1); q_st=reshape(q,I*J,1); if flagN == 0 || flagQ == 0 intN=(gg.*n_st); N1 = sum(intN); % Update aggregate labor intQ=(gg.*q_st); Q1 = sum(intQ); % Update aggregate quantity end
if abs(N-N1)<critN flagN=1; end if abs(Q-Q1)<critQ flagQ=1; end
N = relax*N +(1-relax)*N1; Q = relax*Q +(1-relax)*Q1;
p = pfun(Q); w = wfun(p,zmean,N);
MM=sum(gg(mlocX))/sum(gg);
if flagN==1 && flagQ==1 break end end toc; disp(strcat('exit/entry rate: ',num2str(100*MM),'%')); XR(cfi,cei)=MM; TN(cfi,cei)=N; TQ(cfi,cei)=Q; Tw(cfi,cei)=w; Tp(cfi,cei)=p; endend
%% GRAPHSif length(cfAll)>1 && length(ceAll)>1 figure surf(ceAll,cfAll(1:17),XR(1:17,:)) view([40 50]) ylabel('fixed cost ($c_f$)','FontSize',12,'interpreter','latex') xlabel('entry cost ($c_e$)','FontSize',12,'interpreter','latex') zlabel('exit rate ($\lambda$)','FontSize',12,'interpreter','latex')endTT = (V<-c_en);VV11 = V;VV11(TT) = -c_en;if length(cfAll)==1 && length(ceAll)==1 figure surf(a,z,g') view([45 20]) xlabel('age ($a$)','FontSize',12,'interpreter','latex') ylabel('shocks ($z$)','FontSize',12,'interpreter','latex') zlabel('distribution $f(a,z)$','FontSize',12,'interpreter','latex') figure surf(a,z,VV11') view([45 20]) xlabel('age ($a$)','FontSize',12,'interpreter','latex') ylabel('shocks ($z$)','FontSize',12,'interpreter','latex') zlabel('value $v(a,z)$','FontSize',12,'interpreter','latex')end
% figure;% subplot(2,2,1:2); surf(z,a,pi1); title('profits');% subplot(2,2,3); surf(z,a,n); title('labor');% subplot(2,2,4); surf(z,a,q); title('quantity');