Problem Set 4 questões economia

Prévia do material em texto

Problem Set 4
Jose M Chizzotti
Problem 1 item (a)
% Below follows the codeG linked in the problem set explained
% in detail. In order to differenciate my comments from those
% of the authors of the code I put everything I wrote in CAPS
% and between " ". I explicitly add a paramter for TFP in the
% code where necessary (since this parameter was normalized
% to 1 and so ignored) , in order to investigate the results
% of a increse in this parameter, and I also change the
% initial guess for K accordingly. I assume "hick's neutral"
% technical change, that is Yt = At*F(Kt,Nt), where A is the
% TFP and F(K,N) is the production function.
%SIMULATES AN AIYAGARI ECONOMY WITH IDIOSYNCRATIC BROWNIAN MOTION
%GALO NUNO, based on codes from BENJAMIN MOLL
%Optimized for speed by SeHyoun Ahn
%The algorithm is based on a relaxation scheme to find K. The value
 function is
%found every round by solving the HJB equation through an upwind
 finite
%differences scheme. The distribution is found by also solving using
 finite
%differences the Fokker-Planck (Kolmogorov forward) equation
clear all; close all; format long;
tic;
%--------------------------------------------------
%PARAMETERS
ga = 2; % CRRA utility with parameter gamma
alpha = 0.35; % Production function F = K^alpha * L^(1-alpha)
delta = 0.1; % Capital depreciation
TFP = 1; % "TOTAL FACTOR PRODUCTIVITY (ADDED)"
%--------------------------------------------------
% "HERE THE AUTHORS SET THE PARAMETERS FOR THE DIFFUSION PROCESS
% 'dz_t = mu(z_t)*dt + sigma(z_t)*dW_t', A CONTINUOUS TIME
% GAUSS-MARKOV PROCESS THAT CAN BE SEEN AS THE ANALOUGUE OF THE
% DISCRETE TIME AR(1) PROCESS (ALSO CALLED "ORNSTEIN-UHLENBECK
% PROCESS"). z_t IS THE PRODICTIVITY OF THE WORKER AT PERIOD t."
zmean = 1.0; % mean O-U process (in levels). This parameter has
 to be adjusted to ensure that the mean of z (truncated gaussian) is
 1.
sig2 = (0.10)^2; % sigma^2 O-U
Corr = exp(-0.3); % persistence -log(Corr) O-U
rho = 0.05; % discount rate
1
Problem Set 4
%------------------------------------------------
% "NOTICE THAT THE GUESS USED BY THEM FOR K (3.8) IS SLIGHTLY ABOVE
% THE STATIONARY STATE LEVEL OF K WITH r=rho, WHICH IS Kss =
% (alpha*TFP/(rho+delda))^(1/(1-alpha))~3.68. THAT'S BECAUSE WITH
% PRECAUTIONATY SAVINGS r TENDS TO BE LOWER THAN rho SO IN
% EQUILIBRIUM THE MARGINAL PRODUCTIVITY OF CAPITAL IS LOWER,
% IMPLING HIGHER K. IN ORDER TO ADAPT THIS GUESS FOR A TFP DIFFERENT
% THAN ONE I ADD (dKss/dTFP)*(TFP - 1) TO IT."
% initial aggregate capital. It is important to guess a value close
% to the solution for the algorithm to converge
delta_K = (1/(1-alpha))*((alpha/(rho+delta))^(alpha/(1-
alpha)))*(alpha/(rho+delta))*(TFP-1);
K = 3.8 + delta_K;
relax = 0.99; % relaxation parameter
J=40; % number of z points
zmin = 0.5; % Range z
zmax = 1.5;
amin = -1; % borrowing constraint
amax = 30; % range a
I=100; % number of a points
%simulation parameters
maxit = 100; %maximum number of iterations din the HJB loop
maxitK = 100; %maximum number of iterations in the K loop
crit = 10^(-6); %criterion HJB loop
critK = 1e-5; %criterion K loop
Delta = 1000; %delta in HJB algorithm
%ORNSTEIN-UHLENBECK IN LEVELS
the = -log(Corr);
Var = sig2/(2*the);
%--------------------------------------------------
%VARIABLES
a = linspace(amin,amax,I)'; %wealth vector
da = (amax-amin)/(I-1);
z = linspace(zmin,zmax,J); % productivity vector
dz = (zmax-zmin)/(J-1);
dz2 = dz^2;
aa = a*ones(1,J);
zz = ones(I,1)*z;
%-------------------------------------------------
% "OBTAIN THE PARAMETERS FROM THE O-U PROCESS THAT WILL ENTER IN
% THE HJB AND THE KF EQUATIONS".
mu = the*(zmean - z); %DRIFT (FROM ITO'S LEMMA)
s2 = sig2.*ones(1,J); %VARIANCE (FROM ITO'S LEMMA)
%-------------------------------------------------
%Finite difference approximation of the partial derivatives
2
Problem Set 4
Vaf = zeros(I,J);
Vab = zeros(I,J);
Vzf = zeros(I,J);
Vzb = zeros(I,J);
Vzz = zeros(I,J);
c = zeros(I,J);
%------------------------------------------------
% "HERE THE AUTHORS CONSTRUCT WHAT IS EQUIVALENT TO THE MATRIX B
% IN TIAGO'S LECTURE FOR THE AIYAGARI MODEL WITH 2 STATES (HERE
% THEY ARE APPROXIMTING THE CONTINUOUS PROCESS USING A GRID WITH 40
% STATES). THIS IS AN AUXILIARY MATRIX THAT WILL ENTER IN THE
% MATRICES USED FOR THE LINEAR OPERATIONS IN THE ITERATIONS TO FIND
% THE HJB AND IN THE LINEAR OPERATION TO FIND THE STATIONARY
% DISTRIBUTION OF STATES USING THE KF EQUATIONS. IN TIAGO'S LECTURE
% THESE AUXILIARY MATRICES ARE P_n = S_n*D_n + B, WHERE n IS THE
% INDEX OF THE ITERATION."
%CONSTRUCT MATRIX Aswitch SUMMARIZING EVOLUTION OF z
yy = - s2/dz2 - mu/dz;
chi = s2/(2*dz2);
zeta = mu/dz + s2/(2*dz2);
%This will be the upperdiagonal of the matrix Aswitch
updiag=zeros(I,1); %This is necessary because of the peculiar way
 spdiags is defined.
for j=1:J
 updiag=[updiag;repmat(zeta(j),I,1)];
end
%This will be the center diagonal of the matrix Aswitch
centdiag=repmat(chi(1)+yy(1),I,1);
for j=2:J-1
 centdiag=[centdiag;repmat(yy(j),I,1)];
end
centdiag=[centdiag;repmat(yy(J)+zeta(J),I,1)];
%This will be the lower diagonal of the matrix Aswitch
lowdiag=repmat(chi(2),I,1);
for j=3:J
 lowdiag=[lowdiag;repmat(chi(j),I,1)];
end
% "OBS: SINCE THE MATRIX Aswitch (AND OTHERS BELOW) CONTAINS A LOT
% OF ZEROS, THEY CREATE IT AS A SPARSE BAND MATRIX (MATLAB ONLY
% KEEPS RECORDS OF NON-ZERO ELEMENTS) TO IMRPOVE EFFICIENCY OF THE
% CODE."
%Add up the upper, center, and lower diagonal into a sparse matrix
Aswitch=spdiags(centdiag,0,I*J,I*J)+spdiags(lowdiag,-
I,I*J,I*J)+spdiags(updiag,I,I*J,I*J);
%----------------------------------------------------
3
Problem Set 4
% "USING THEIR INITIAL GUESS FOR K THEY SET INITIAL r TO EQUAL THE
% MARGINAL PRODUCTIVITY OF CAPITAL - DEPRECIATION AND w TO EQUAL
% MARGINAL PRODUCTIVITY OF LABOR IMPLIED BY K. THE INITIAL GUESS
% FOR THE VALUE FUNCTION IS THAT IMPLIED BY ZERO SAVINGS (STEADY
% STATE) AND ZERO PROBABILITY OF CHANGING PRODUCTIVITY STATE IN
% EACH POINT OF THE STATES GRID. I ADD THE TFP (THAT WAS
% IMPLICITLY = 1) IN THE EQUATIONS FOR WAGE AND INTEREST RATE"
%INITIAL GUESS
r = TFP*(alpha * K^(alpha-1)) -delta; %interest rates "TFP ADDED"
w = TFP*((1-alpha) * K^(alpha)); %wages "TFP ADDED"
v0 = (w*zz + r.*aa).^(1-ga)/(1-ga)/rho;
v = v0;
dist = zeros(1,maxit);
%-----------------------------------------------------
% "STARTS ITERATION TO FIND PRICES "r" AND "w" AND AGGREGATE CAPITAL
% USING RELAXATION ALGORITHM".
%MAIN LOOP
for iter=1:maxitK
%disp('Main loop iteration')
%disp(iter)
%-----------------------------------------------------
% "STARTS ITERATION TO FIND THE HJB EQUATION USING AN UPWIND FINITE
% DIFFERENCE SCHEME. THAT IS, FOR POINTS ON THE STATE GRID WHERE THE
% DRIFT IN THE STATE VARIABLE a IS POSITVE(NEGATIVE), (OPTIMAL SAVINGS
% IMPLIED BY c = u'^-1(V_a(a,z)) IS POSITIVE (NEGATIVE)), THE
% DERIVATIVE OF THE VALUE FUNCTION WITH RESPECT TO a IS APPROXIMATED
% USING THE FORWARD (BACKWARD) DIFFERENCE."
 % HAMILTON-JACOBI-BELLMAN EQUATION %
 for n=1:maxit
 V = v;
 % forward difference
 Vaf(1:I-1,:) = (V(2:I,:)-V(1:I-1,:))/da;
 Vaf(I,:) = (w*z + r.*amax).^(-ga); %will never be used, but
 impose state constraint a Vaf; %indicator whether value
 function is concave (problems arise if this is not the case)
 %consumption and savings with forward difference
 cf = Vaf.^(-1/ga);
 sf = w*zz + r.*aa - cf;
 %consumption and savings with backward difference
4
Problem Set 4
 cb = Vab.^(-1/ga);
 sb = w*zz + r.*aa - cb;
 %consumption and derivative of value function at steady state
 c0 = w*zz + r.*aa;
 Va0 = c0.^(-ga);
 % "HERE THEY EMPLOY THE UPWINDSCHEME TO APPROXIMATE dV/da AND
 % GET OPTIMAL c IMPLIED BU V_n, WHICH WILL BE USED TO OBTAIN
 % THE UPDATE V_n+1"
 % dV_upwind makes a choice of forward or backward differences
 based on
 % the sign of the drift
 If = sf > 0; %positive drift --> forward difference
 Ib = sb backward difference
 I0 = (1-If-Ib); %at steady state
 %make sure backward difference is used at amax
 % Ib(I,:) = 1; If(I,:) = 0;
 %STATE CONSTRAINT at amin: USE BOUNDARY CONDITION UNLESS sf >
 0:
 %already taken care of automatically
 Va_Upwind = Vaf.*If + Vab.*Ib + Va0.*I0; %important to include
 third term
 c = Va_Upwind.^(-1/ga);
 u = c.^(1-ga)/(1-ga);
 %--------------------------------------------
 % "MATRIX AA IS EQUIVALENT TO MATRIX S_n*D_n IN TIAGO'S
 % LECTURE AND MATRIX A IS EQUVALENT TO P_n."
 %CONSTRUCT MATRIX A
 X = - min(sb,0)/da;
 Y = - max(sf,0)/da + min(sb,0)/da;
 Z = max(sf,0)/da;
 updiag=[0]; %This is needed because of the peculiarity of
 spdiags.
 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);
5
Problem Set 4
 % "EQUIVALENT TO P_n = S_n*D_n + B."
 A = AA + Aswitch;
 %-----------------------------------------
 % "HERE THEY USE THE "IMPLICIT METHOD" FOR THE VALUE FUNTION
 % ITERATION, IN WHICH V_n+1 IS OBTAINED BY PREMULTIPLICATING
 % u(c) + (1/Delta)*V_n BY ((1/Delta + rho)*I - A_n)^(-1),
 % WHERE c IS OPTIMAL CONSUMPTION AND A_n IS THE MATRIX
 % CONSTRUCTED ABOVE, BOTH A FUNCTION OF V_n. IN HIS LECTURE
 % MOLL ARGUES THIS METHOD IS PREFERABLE OVER THE EXPLICIT
 % METHOD, WHICH JUST SUBSTITUTE THE OPTIMAL CONSUMPTION
 % IMPLIED BY V_n IN THE HJB, CALCULATES THE DIFFERENCE
 % BETWEEN THIS VALUE AND rho*V_n AND GET THE UPDATE AS
 % V_n+1 = V_n + Delta*DIFFERENCE. THIS IS BEACUASE THE
 % IMPLICIT METHOD IS GENERALY MORE EFFICIENT AND STABLE.
 % IN BOTH METHODS, Delta IS A PARAMETER THAT GIVES THE
 % SPEED OF THE UPDATE, AND GENERALLY, IT HAS TO BE A SMALL
 % VALUE SINCE THE CONTRACTION MAPPING THEOREM DOES NOT
 % WORK IN CONTINUOUS TIME, SO CONVERGENCE IS NOT GUARANTEED.
 % (0BS: SINCE u AND V_n ENTER IN THE LINEAR OPARATION AS
 % VECTORS, AND THE AUTHOS CRONSTRUCTED THEM AS MATRICES,
 % THESE OBJECTS NEED TO BE STACKED ACCORDINGLY)."
 B = (1/Delta + rho)*speye(I*J) - A;
 u_stacked = reshape(u,I*J,1);
 V_stacked = reshape(V,I*J,1);
 b = u_stacked + V_stacked/Delta;
 % "HERE "B\b" STANDS FOR THE MATRIX OPERATION (B^(-1))*b"
 V_stacked = B\b; %SOLVE SYSTEM OF EQUATIONS
 V = reshape(V_stacked,I,J);
 Vchange = V - v;
 v = V;
 dist(n) = max(max(abs(Vchange)));
 if dist(n)rate: 0.0486 
Equilibrium aggregate capital: 3.7365 
Equilibrium GDP: 1.5862 
Equilibrium average consumption: 1.2126 
Equilibrium capital to output ratio: 2.3556 
Equilibrium mean of wealth distribution: 3.7364 
Equilibrium standard deviation of wealth distribution: 3.5571 
Equilibrium median of wealth distribution: 3.0707 
Equilibrium mean of productivity distribution: 1.0001 
Equilibrium standard deviation of productivity distribution: 0.1289 
9
Problem Set 4
10
Problem Set 4
Problem 1 items (b) and (c)
% In order to performe the comparative statics exercise I put the
% script above in a .m file as a funtion that takes as input TFP and
% returns all the statistics analized in item (a) plus the graphs.
% Since I kept the fprintfs and the graphs in the function's code I
% only have to call it. (To clarify why I did not make the code
% above already as a function, that is because functions have to
% go in the end of the script when publishing with Matlab)
[K2,Delta_K2,r2]=solve_aiy(1.1);
fprintf('\n')
fprintf('Change in K approximated using derivative %.4f \n',Delta_K2)
fprintf('Actual change in K %.4f \n',K2-K)
fprintf('Change in interest rate %d \n',r2-r)
% First of all, since stationaty state capital depends positively
% on Total Factors Productivity we would expect this to raise by
% around (dK_ss/dTFP)*(0.1) = 0.5694 (this is given by Delta_K),
% which is indeed the case, since K2 - K = 0.5898. Since wages
% depend positively on TFP and capital, wages should increase
% and the same for output (since labor suply is inelastic in this
% model), the resuts are in line with this reasoning. What is
% interest to notice is that the interest rate stay the same
% (r2-r) = 4e-6, which given that we are approximating solutions
% with margins around this magnitude, is not significantly
% different from zero). That is, the positive effect of the higher
% TFP couterbalances the negative effect from higher capital. The
% intuition is that the households are still subject to the same
% shocks as before, so assuming a fixed interst rate, they raise
% their precautionary savings in proportion to the increase of their
% income, and since non-precautionaty savings (savings in the absence
% of uncertainty) also increase in proportion to income, total
% savings, therefore asset holdings in equilibrium, increase in
% proportion to income. Now, we know that for a fixed interest rate
% (and inelastic labor supply), with the specified production
% function, in equilibrium, output increase in the same proportion
% as capital demand with a increase in the TFP, so income also
% increase in proportion to capital demand in equlibrium. Hence,
% capital supply and capital demand increase in the same proportion
% and the interst rate before the TFP shock still clears the market.
% The capital to output stays the same and this was expected by the
% reasoning above. We can see that inequality in wealth
% distribution increases, since there is an increase in the standart
% deviation and in the difference between mean and median of the
% distribution. Finally, average consumption increases in the same
% proportion as capital and output.
New equilibrium wage rate: 1.1938 
New equilibrium interest rate: 0.0486 
New equilibrium aggregate capital: 4.3263 
New equilibrium GDP: 1.8367 
New equilibrium average consumption: 1.4041 
11
Problem Set 4
New equilibrium capital to output ratio: 2.3555 
New equilibrium mean of wealth distribution: 4.3263 
New equilibrium standard deviation of wealth distribution: 4.0078 
New equilibrium median of wealth distribution: 3.3838 
Equilibrium mean of productivity distribution: 1.0001 
Equilibrium standard deviation of productivity distribution: 0.1289 
Change in K approximated using derivative 0.5694 
Actual change in K 0.5899 
Change in interest rate 4.339804e-06 
12
Problem Set 4
Problem 1 item (d)
% Ramsey-Cass-Koopmans Model (I am assuming that the item only asks us
% to reproduce the comparative statics from item (c), therefore I will
% only compute the equilibirum values for the variables).
alp = 1; % avarage labor productivity, the same as in aiyagari model
 % above (zmean=1);
r_rck = rho;
K_rck = (TFP*alpha/(rho + delta))^(1/(1-alpha));
w_rck = TFP*((1-alpha) * K_rck^(alpha));
gdp_rck = TFP*K_rck^alpha*alp^(1-alpha);
avg_c_rck = gdp_rck - delta*K_rck;
ky_ratio_rck = K_rck/gdp_rck;
mean_a_rck = K_rck;
sd_a_rck = 0;
fprintf('Equilibrium wage rate: %.4f \n',w_rck);
fprintf('Equilibrium interest rate: %.4f \n',r_rck);
fprintf('Equilibrium aggregate capital: %.4f \n',K_rck);
fprintf('Equilibrium GDP: %.4f \n',gdp_rck);
fprintf('Equilibrium average consumption: %.4f \n',avg_c_rck)
fprintf('Equilibrium capital to output ratio: %.4f \n',ky_ratio_rck);
fprintf('Equilibrium mean of wealth distribution: %.4f
 \n',mean_a_rck);
13
Problem Set 4
fprintf('Equilibrium standard deviation of wealth distribution: %.4f
 \n',sd_a_rck);
fprintf('\n')
solve_rck(1.1)
% The signs and magnitutes of the impact from the increase in TFP in
% the variables of the model are the same as in the Aiyagari model,
% except for the sd of wealth distribution which now is zero since
% all households are identical and therefore have the same wealth in
% equilibrium. What is intersting to observe are the differences
% between the variables in this model's equilibrium compared
% to the same variables in the equilibrium from Aiyagari's model.
% In Aiyagari's model, the introduction of uncertainty and incomplete
% markets cause the households to engage in precautionaty savings in
% order to partially insure themselfs against the idiossincratic
% shock to which they are subjected. This means that for given
% parameters and prices, there will be more saving by the households
% and therefore, comapared to the Ramsey model, aggregate capital
% will be higher and the interest rate will be lower in order to
% clear the assets' market. That implies that output and wages will
% be higher and, since the economy is never dynamically inefficient
% in the Ramsey model (capital level is below golden rule),
% consumption will also increase. Finally, the lower interst rate
% implies a higher capital to output ratio in equilibrium, since the
% marginal productivity of capital will be lower. So we can see that
% the differences in the variables computed above and those from
% items (a-c) are in line with the theory.
Equilibrium wage rate: 1.0258 
Equilibrium interest rate: 0.0500 
Equilibrium aggregate capital: 3.6823 
Equilibrium GDP: 1.5781 
Equilibrium average consumption: 1.2099 
Equilibrium capital to output ratio: 2.3333 
Equilibrium mean of wealth distribution: 3.6823 
Equilibrium standard deviation of wealth distribution: 0.0000 
New equibrium wage rate: 1.1878 
New equibrium interest rate: 0.0500 
New equibrium aggregate capital: 4.2638 
New equibrium GDP: 1.8274 
New equilibrium average consumption: 1.4010 
New equibrium capital to output ratio: 2.3333 
New equibrium mean of wealth distribution: 4.2638 
New equibrium standard deviation of wealth distribution: 0.0000 
Published with MATLAB® R2020b
14
Problem 2 item (a)
% Set parameters of AR(1) process
rho=0.98; %autoregressive coefficient
sigma2_z = 0.621; %variance of AR(1)
sigma_e = sqrt(sigma2_z*(1-rho^2)); %standard deviation of white
 %noise
Nz = 7; %Number of states in approximation using markov process
% Calculates the states vector, transition matrix and invariant
% long run distribution
[Pi,ln_z,invdist_z]=markovappr(rho,sigma_e,3,Nz);
z=exp(ln_z);
% Simulates markov chain and plots the results
mkv_chain = markovsimul(Pi,z,1000,4,12);
plot(mkv_chain)
Published with MATLAB® R2020b
1
Problem 2 item (c) and (d)
% Paramters
beta = 0.94; % discount factor
mu=2; % CRRA coeficient
w=1; % wage (c is assumed to be normalized so that w=1)
rh = 0.05; % interest rate of high return asset
rl = 0.01; % interest rate of low return asset
% Grid for asssets
a_max=100;% maximum assets
a_min=0; % minimum assets (no borrowing)
Na=2000; % Number of points in the grid
a=linspace(a_min,a_max,Na)';
psi_max = 4; % Initial maximum psi ==> a_bar ~ 100
psi_min = 0.001; % Initial minimum psi
err = 1;
tol = 0.0001;
iter = 0;
% Start iteration to find psi using bissection algorithm
while err > tol
iter = iter+1;
psi = (psi_max+psi_min)/2;
% Quantity of investment that makes households indifferent between
% assets
a_bar = ((1+rh)/(rh-rl))*psi;
% Calculates consumption and utility grid
c = zeros(Na,Na,Nz);
u = zeros(Na,Na,Nz);
for ii = 1:Na
 for jj = 1:Na
 for kk = 1:Nz
 a1 = a(ii); % assets tomorrow
 a0 = a(jj); % asets today
 z0 = z(kk); % productivity
 % Calculates consumption taking into consideration if
 % the agent chose the high or the low return asset last
 % period
 if a0=a_bar
 c(ii,jj,kk)= (1+rh)*(a0-psi) +w*z0 -a1;
 end
1
 % Drop points where consumption is negative from
 % consideration
 if c(ii,jj,kk) vftol
 vfit=vfit+1;
 V1 = TV1;
 V2 = TV2;
 V3 = TV3;
 V4 = TV4;
 V5 = TV5;
 V6 = TV6;
 V7 = TV7;
 % expected value function
 aux1 =
 Pi(1,1)*V1+Pi(1,2)*V2+Pi(1,3)*V3+Pi(1,4)*V4+Pi(1,5)*V5+Pi(1,6)*V6+Pi(1,7)*V7;
 aux2 =
 Pi(2,1)*V1+Pi(2,2)*V2+Pi(2,3)*V3+Pi(2,4)*V4+Pi(2,5)*V5+Pi(2,6)*V6+Pi(2,7)*V7;
 aux3 =
 Pi(3,1)*V1+Pi(3,2)*V2+Pi(3,3)*V3+Pi(3,4)*V4+Pi(3,5)*V5+Pi(3,6)*V6+Pi(3,7)*V7;
 aux4 =
 Pi(4,1)*V1+Pi(4,2)*V2+Pi(4,3)*V3+Pi(4,4)*V4+Pi(4,5)*V5+Pi(4,6)*V6+Pi(4,7)*V7;
 aux5 =
 Pi(5,1)*V1+Pi(5,2)*V2+Pi(5,3)*V3+Pi(5,4)*V4+Pi(5,5)*V5+Pi(5,6)*V6+Pi(5,7)*V7;
 aux6 =
 Pi(6,1)*V1+Pi(6,2)*V2+Pi(6,3)*V3+Pi(6,4)*V4+Pi(6,5)*V5+Pi(6,6)*V6+Pi(6,7)*V7;
2
 aux7 =
 Pi(7,1)*V1+Pi(7,2)*V2+Pi(7,3)*V3+Pi(7,4)*V4+Pi(7,5)*V5+Pi(7,6)*V6+Pi(7,7)*V7;
 TV1 = (max(u(:,:,1) + beta*aux1*o))';
 TV2 = (max(u(:,:,2) + beta*aux2*o))';
 TV3 = (max(u(:,:,3) + beta*aux3*o))';
 TV4 = (max(u(:,:,4) + beta*aux4*o))';
 TV5 = (max(u(:,:,5) + beta*aux5*o))';
 TV6 = (max(u(:,:,6) + beta*aux6*o))';
 TV7 = (max(u(:,:,7) + beta*aux7*o))';
 vferr=max([max(abs(TV1-V1)),max(abs(TV2-V2)),max(abs(TV3-V3)),...
 max(abs(TV4-V4)),max(abs(TV5-V5)),max(abs(TV6-
V6)),max(abs(TV7-V7))]);
end
% Get the indices of policy funtion for assets
[TV1,i1] = max(u(:,:,1) + beta*aux1*o);
[TV2,i2] = max(u(:,:,2) + beta*aux2*o);
[TV3,i3] = max(u(:,:,3) + beta*aux3*o);
[TV4,i4] = max(u(:,:,4) + beta*aux4*o);
[TV5,i5] = max(u(:,:,5) + beta*aux5*o);
[TV6,i6] = max(u(:,:,6) + beta*aux6*o);
[TV7,i7] = max(u(:,:,7) + beta*aux7*o);
% policy function for assets matrix
policy_a = [a(i1), a(i2), a(i3), a(i4), a(i5), a(i6),a(i7)];
% Constructs policy function for consumption based on policy function
% for assets
policy_c = zeros(Na,Nz);
for i = 1:Na
 for j = 1:Nz
 a0 = a(i);
 a1 = policy_a(i,j);
 z0 = z(j);
 if a0=a_bar
 policy_c(i,j)= (1+rh)*(a0-psi) +w*z0 -a1;
 end
 end
end
% Calculates endogenous markov chain in order to start simulation
% below with wealth and productitvity distributions for households
% that reflect stationary state distribution
% Constructs matrices with index of policy funtion for each
% state
proba1=zeros(Na,Na);
3
proba2=zeros(Na,Na);
proba3=zeros(Na,Na);
proba4=zeros(Na,Na);
proba5=zeros(Na,Na);
proba6=zeros(Na,Na);
proba7=zeros(Na,Na);
for i=1:Na
 for j =1:Na
 if i1(i)==j
 proba1(i,j)=1;
 else
 proba1(i,j)=0;
 end
 if i2(i)==j
 proba2(i,j)=1;
 else
 proba2(i,j)=0;
 end
 if i3(i)==j
 proba3(i,j)=1;
 else
 proba3(i,j)=0;
 end
 if i4(i)==j
 proba4(i,j)=1;
 else
 proba4(i,j)=0;
 end
 if i5(i)==j
 proba5(i,j)=1;
 else
 proba5(i,j)=0;
 end
 if i6(i)==j
 proba6(i,j)=1;
 else
 proba6(i,j)=0;
 end
 if i7(i)==j
 proba7(i,j)=1;
 else
 proba7(i,j)=0;
 end
 end
end
% Calculates Matrix of transition for assets and productivity
proba = zeros(Na*Nz,Na*Nz);
for i=1:Na
 for j=1:Na
 proba(7*i-6,7*j-6)=Pi(1,1)*proba1(i,j);
 proba(7*i-6,7*j-5)=Pi(1,2)*proba1(i,j);
4
 proba(7*i-6,7*j-4)=Pi(1,3)*proba1(i,j);
 proba(7*i-6,7*j-3)=Pi(1,4)*proba1(i,j);
 proba(7*i-6,7*j-2)=Pi(1,5)*proba1(i,j);
 proba(7*i-6,7*j-1)=Pi(1,6)*proba1(i,j);
 proba(7*i-6,7*j) =Pi(1,7)*proba1(i,j);
 proba(7*i-5,7*j-6)=Pi(2,1)*proba2(i,j);
 proba(7*i-5,7*j-5)=Pi(2,2)*proba2(i,j);
 proba(7*i-5,7*j-4)=Pi(2,3)*proba2(i,j);
 proba(7*i-5,7*j-3)=Pi(2,4)*proba2(i,j);
 proba(7*i-5,7*j-2)=Pi(2,5)*proba2(i,j);
 proba(7*i-5,7*j-1)=Pi(2,6)*proba2(i,j);
 proba(7*i-5,7*j) =Pi(2,7)*proba2(i,j);
 proba(7*i-4,7*j-6)=Pi(3,1)*proba3(i,j);
 proba(7*i-4,7*j-5)=Pi(3,2)*proba3(i,j);
 proba(7*i-4,7*j-4)=Pi(3,3)*proba3(i,j);
 proba(7*i-4,7*j-3)=Pi(3,4)*proba3(i,j);
 proba(7*i-4,7*j-2)=Pi(3,5)*proba3(i,j);
 proba(7*i-4,7*j-1)=Pi(3,6)*proba3(i,j);
 proba(7*i-4,7*j) =Pi(3,7)*proba3(i,j);
 proba(7*i-3,7*j-6)=Pi(4,1)*proba4(i,j);
 proba(7*i-3,7*j-5)=Pi(4,2)*proba4(i,j);
 proba(7*i-3,7*j-4)=Pi(4,3)*proba4(i,j);
 proba(7*i-3,7*j-3)=Pi(4,4)*proba4(i,j);
 proba(7*i-3,7*j-2)=Pi(4,5)*proba4(i,j);
 proba(7*i-3,7*j-1)=Pi(4,6)*proba4(i,j);
 proba(7*i-3,7*j) =Pi(4,7)*proba4(i,j);
 proba(7*i-2,7*j-6)=Pi(5,1)*proba5(i,j);
 proba(7*i-2,7*j-5)=Pi(5,2)*proba5(i,j);
 proba(7*i-2,7*j-4)=Pi(5,3)*proba5(i,j);
 proba(7*i-2,7*j-3)=Pi(5,4)*proba5(i,j);
 proba(7*i-2,7*j-2)=Pi(5,5)*proba5(i,j);
 proba(7*i-2,7*j-1)=Pi(5,6)*proba5(i,j);
 proba(7*i-2,7*j) =Pi(5,7)*proba5(i,j);
 proba(7*i-1,7*j-6)=Pi(6,1)*proba6(i,j);
 proba(7*i-1,7*j-5)=Pi(6,2)*proba6(i,j);
 proba(7*i-1,7*j-4)=Pi(6,3)*proba6(i,j);
 proba(7*i-1,7*j-3)=Pi(6,4)*proba6(i,j);
 proba(7*i-1,7*j-2)=Pi(6,5)*proba6(i,j);
 proba(7*i-1,7*j-1)=Pi(6,6)*proba6(i,j);
 proba(7*i-1,7*j) =Pi(6,7)*proba6(i,j);
 proba(7*i,7*j-6)=Pi(7,1)*proba7(i,j);
 proba(7*i,7*j-5)=Pi(7,2)*proba7(i,j);
 proba(7*i,7*j-4)=Pi(7,3)*proba7(i,j);
 proba(7*i,7*j-3)=Pi(7,4)*proba7(i,j);
 proba(7*i,7*j-2)=Pi(7,5)*proba7(i,j);
 proba(7*i,7*j-1)=Pi(7,6)*proba7(i,j);
 proba(7*i,7*j) =Pi(7,7)*proba7(i,j);
5
 end
end
% Iterates using the endogenous transition matrix in order to obtain
% long run distribution of assets and productivity
maxit=1000;
I0 = ones(1,Na*Nz)*1/(Na*Nz);
for i=1:maxit
 I=I0*proba;
 if norm(I-I0)=maxit
 fprintf('WARNING: Maximum number of % iterations reached',maxit)
end
% Calculates the stationary marginal distribution of assets
margdist_a = zeros(Na,1);
for i=1:Na
 margdist_a(i) = sum(I(7*i-6:7*i));
end
cdf_a = zeros(Na,1);
cdf_a(1) = margdist_a(1);
for i=2:Na
 cdf_a(i) = cdf_a(i-1) + margdist_a(i);
end
cdf_z = zeros(Nz,1);
cdf_z(1) = invdist_z(1);
for i=2:Nz
 cdf_z(i) = cdf_z(i-1)+ invdist_z(i);
end
N = 1000; % number of households simulated
T=50; % number of periods for the simulation
wealth = zeros(N,20);
income = zeros(N,20);
for i=1:N
 rng(i); % set seed, so that the same random numbers are used for
 % in each iteration for psi.
 % Generates a random draw from the stationary distribution of a
 rd1 = rand;
 aux_a = find((cdf_a - rd1)>0);
 ind = aux_a(1);
 a0=a(ind);
6
 % Generates a random dram from the stationary distribution of z
 rd2 = rand;
 aux_z = find((cdf_z - rd2)>0);
 ind2 = aux_z(1);
 z0 = z(ind2);
 % Simulates a markov chain for the household productivity
 % evolution in time
 [chain_z,state_z] = markovsimul(Pi,z,T,ind2,i);
 path_a = zeros(T,1); % path of assets
 path_a(1)=a0; % initial assets
 path_inc = zeros(T,1); % path of income
 % Calculates initial income based on a0
 if a0 a_bar;
quant_highr = sum(aux_m(:));
prop_highr = quant_highr/(N*20);
% implements bissection algorithm
if prop_highr > 0.1 %if to many households are investing in the high
 psi_min=psi; % return asset, increase cost psi
else % otherwise, if too little are investing
 psi_max=psi; % decrease cost psi
end
err = psi_max-psi_min; % iterate until bissection has converged to
7
 % somewhere and check if the value obtained
 % is indeed the root.
if prop_highr == 0.1
 break
end
end
% Report the results
fprintf('Value function iterations: %i \n',vfit)
fprintf('itereations for psi: %i \n',iter)
fprintf('Calibrated psi: %.4f \n',psi)
fprintf('Proportion of households saving on high asset: %.4f%%
 \n',prop_highr*100)
fprintf('Quantity of assets that make households indiffent between
 assets: %.4f \n',a_bar);
fprintf('\n')
% Plots policy functions for highest productivity state
figure(1)
plot(a,policy_a(:,7));
title('policy function for investment (highest prod)');
xlabel('assets today');
ylabel('assets tomorrow');
figure(2)
plot(a,policy_c(:,7))
title('policy function for consumption (highest prod)');
xlabel('assets');
ylabel('consumption');
Value function iterations: 169 
itereations for psi: 13 
Calibrated psi: 0.5355 
Proportion of households saving on high asset: 10.0000% 
Quantity of assets that make households indiffent between assets:
 14.0578 
8
9
W1 = wealth(:);
% Estimates stationary pdf of wealth using data from simulation
pdf_w = zeros(Na,1);
for i=1:Na
 aux_pdfw = W1==a(i);
 pdf_w(i) = sum(aux_pdfw)/(N*20);
end
% Estimates stationary cdf of wealth using data from simulation
cdf_w = zeros(Na,1);
cdf_w(1) = pdf_w(1);
for i=2:Na
 cdf_w(i) = cdf_w(i-1) + pdf_w(i);
end
% Repeat the process for the distribution of wealth conditional
% on wealth being strictly positive
aux_W2 = find(W1>0);
W2 = W1(aux_W2);
pdf_w2 = zeros(Na,1);
for i=1:Na
 aux_pdfw2 = W2==a(i);
 pdf_w2(i) = sum(aux_pdfw2)/length(W2);
end
cdf_w2 = zeros(Na,1);
cdf_w2(1) = pdf_w2(1);
for i=2:Na
 cdf_w2(i) = cdf_w2(i-1) + pdf_w2(i);
end
% Estimates income shares
total_w = sum(W1);
aux_ineq = sort(W1);
aux_ineq01 = aux_ineq(end-19:end);
share01 = sum(aux_ineq01)/total_w;
aux_ineq1 = aux_ineq(end-199:end);
share1 = sum(aux_ineq1)/total_w;
aux_ineq5 = aux_ineq(end-999:end);
share5 = sum(aux_ineq5)/total_w;
aux_ineq10 = aux_ineq(end-1999:end);
share10 = sum(aux_ineq10)/total_w;
aux_ineq50 = aux_ineq(end-9999:end);
share50 = sum(aux_ineq50)/total_w;
% Gini coefficient for wealth
aux_gini = zeros(N*20,N*20);
10
for i=1:N*20
 for j=1:N*20
 aux_gini(i,j) = abs(W1(i)-W1(j));
 end
end
aux_gini2 = aux_gini(:);
mean_w = mean(W1);
gini = sum(aux_gini2)/(2*(N*20)^2*mean_w);
sd_w = std(W1); % standard deviation of wealth dist
median_w = median(W1); % median of wealth dist
% Report estimated statistics for the wealth distribution
fprintf('Wealth distributon mean: %.4f \n',mean_w)
fprintf('Wealth distributon standard deviation: %.4f \n',sd_w)
fprintf('Wealth distributon median: %.4f \n',median_w)
fprintf('Top 0.1%% share of wealth: %.4f%% \n',100*share01)
fprintf('Top 1%% share of wealth: %.4f%% \n',100*share1)
fprintf('Top 5%% share of wealth: %.4f%% \n',100*share5)
fprintf('Top 10%% share of wealth: %.4f%% \n',100*share10)
fprintf('Top 50%% share of wealth: %.4f%% \n',100*share50)
fprintf('Maximum wealth: %.4f \n',aux_ineq(end))
fprintf('Minimum wealth: %.4f \n',aux_ineq(1))
fprintf('GINI coefficient for wealth: %.4f \n',gini)
fprintf('\n')
% Plots figures for wealth distribution
figure(1)
plot(a,pdf_w);
title('Wealth distribution PDF');
xlabel('Wealth');
figure(2)
plot(a,cdf_w); ylim([0 1]);
title('Wealth distribution CDF');
xlabel('Wealth');
figure(3)
histogram(W1);
title('Wealth distribution histogram');
xlabel('Wealth');
% Plots figures for wealth distribution conditional on strictly
% positive wealth
figure(4)
plot(a,pdf_w2);
title('Wealth distribution PDF (for postive wealth)');
xlabel('Wealth');
figure(5)
plot(a,cdf_w2); ylim([0 1]);
title('Wealth distribution CDF (for postive wealth)');
xlabel('Wealth');
11
figure(6)
histogram(W2);
title('Wealth distribution histogram (for postive wealth)');
xlabel('Wealth');
Wealth distributon mean: 4.6031 
Wealth distributon standard deviation: 14.7152 
Wealth distributon median: 0.0000 
Top 0.1% share of wealth: 1.7885% 
Top 1% share of wealth: 16.9235% 
Top 5% share of wealth: 65.2758% 
Top 10% share of wealth: 98.1058% 
Top 50% share of wealth: 100.0000% 
Maximum wealth: 84.8424 
Minimum wealth: 0.0000 
GINI coefficient for wealth: 0.9204 
12
13
14
% Reapeat all the steps done above for income
INC1 = income(:);
INCC = unique(INC1);
pdf_inc = zeros(length(INCC),1);
for i=1:length(INCC)
 aux_pdfinc = INC1==INCC(i);
 pdf_inc(i) = sum(aux_pdfinc)/(N*20);
end
cdf_inc = zeros(length(INCC),1);
cdf_inc(1) = pdf_inc(1);
for i=2:length(INCC)
 cdf_inc(i) = cdf_inc(i-1) + pdf_inc(i);
end
total_inc = sum(INC1);
aux_ineqinc = sort(INC1);
aux_ineqinc01 = aux_ineqinc(end-19:end);
shareinc01 = sum(aux_ineqinc01)/total_inc;
aux_ineqinc1 = aux_ineqinc(end-199:end);
shareinc1 = sum(aux_ineqinc1)/total_inc;
aux_ineqinc5 = aux_ineqinc(end-999:end);
15
shareinc5 = sum(aux_ineqinc5)/total_inc;
aux_ineqinc10 = aux_ineqinc(end-1999:end);
shareinc10 = sum(aux_ineqinc10)/total_inc;
aux_ineqinc50 = aux_ineqinc(end-9999:end);
shareinc50 = sum(aux_ineqinc50)/total_inc;
% Gini coefficient for income
aux_giniinc = zeros(N*20,N*20);
for i=1:N*20
 for j=1:N*20
 aux_giniinc(i,j) = abs(INC1(i)-INC1(j));
 end
end
aux_giniinc2 = aux_giniinc(:);
mean_inc = mean(INC1);
giniinc = sum(aux_giniinc2)/(2*(N*20)^2*mean_inc);
sd_inc = std(INC1);
median_inc = median(INC1);
% Plot the estimated statistics for income distribution
fprintf('Income distributon mean: %.4f \n',mean_inc)
fprintf('Income distributon standard deviation: %.4f \n',sd_inc)
fprintf('Income distributon median: %.4f \n',median_inc)
fprintf('Top 0.1%% share of income: %.4f%% \n',100*shareinc01)
fprintf('Top 1%% share of income: %.4f%% \n',100*shareinc1)
fprintf('Top 5%% share of income: %.4f%% \n',100*shareinc5)
fprintf('Top 10%% share of income: %.4f%% \n',100*shareinc10)fprintf('Top 50%% share of income: %.4f%% \n',100*shareinc50)
fprintf('Maximum income: %.4f \n',aux_ineqinc(end))
fprintf('Minimum income: %.4f \n',aux_ineqinc(1))
fprintf('GINI coefficient for income: %.4f \n',giniinc)
% Plot figures for income distribution
figure(7)
plot(INCC,pdf_inc);
title('Income distribution PDF');
xlabel('Income');
figure(8)
plot(INCC,cdf_inc); ylim([0 1]);
title('Income distribution CDF');
xlabel('Income');
figure(9)
histogram(INC1);
title('Income distribution histogram');
xlabel('Income');
Income distributon mean: 1.9375 
Income distributon standard deviation: 2.1794 
16
Income distributon median: 1.0000 
Top 0.1% share of income: 0.7518% 
Top 1% share of income: 5.8690% 
Top 5% share of income: 23.5013% 
Top 10% share of income: 35.9813% 
Top 50% share of income: 84.2370% 
Maximum income: 14.8499 
Minimum income: 0.0940 
GINI coefficient for income: 0.5127 
17
18
Problem 2 items (e) and (f) part one
% In order to solve for the economy with only one asset, I took the
% script above, set psi_max and psi_min to zero, a_bar to 200 and
% made it as a function on separate .m file that takes as input r,
% sets rh=rl=r, and gives as output the same as above.
solvep2(0.015)
% At first I solved the model for r=0.012, but at interst rates
% this low nobady saves anything and aggregate wealth is zero.
% so I used r=0.015 instead. We would expect inequality in
% wealth and income disribution to decrease now that rich
% households do not have acess to different returns. But this,
% does not heppens. The problem is that with such a low interest
% rate (notice 1/beta -1 ~ 0.63), only those with the highest
% productivity, a very low percentage of the population, engage
% in precautionary savings, so that the inequality of wealth and
% income distribution is worse than before when measured by
% the top percentils shares and the GINI coeficient. Notice
% that the standard deviation of wealth is greatly reduced though,
% as is maximum wealth.
% In order to investigate the effects of giving the poor households
% access to the high return asset, without the fixed cost of
% investment, I ran the same experiment but now with r=0.05.
Value function iterations: 169 
Wealth distributon mean: 0.0104 
Wealth distributon standard deviation: 0.0589 
Wealth distributon median: 0.0000 
Top 0.1% share of wealth: 3.3573% 
Top 1% share of wealth: 33.5731% 
Top 5% share of wealth: 100.0000% 
Top 10% share of wealth: 100.0000% 
Top 50% share of wealth: 100.0000% 
Maximum wealth: 0.3502 
Minimum wealth: 0.0000 
GINI coefficient for wealth: 0.9702 
Income distributon mean: 1.7136 
Income distributon standard deviation: 2.0560 
Income distributon median: 1.0000 
Top 0.1% share of income: 0.6209% 
Top 1% share of income: 6.2089% 
Top 5% share of income: 24.4906% 
Top 10% share of income: 38.6007% 
Top 50% share of income: 83.5716% 
Maximum income: 10.6398 
Minimum income: 0.0940 
GINI coefficient for income: 0.5191 
19
20
21
22
23
Problem 2 items (e) and (f) part two
% Model with one asset, psi=0 and r=0.05
solvep2(0.05)
% Now, not only does avarage wealth increases greatly, without a
% major increase in standard deviation, but also wealth and income
% distribution inequality is greatly reduced, both as measured by
% the top percentils shares and the GINI coefficient. Also, notice
% that the medians for the wealth and income distribution are now
% much closer to its respectives means (in percentage terms),
% another indication of a reduction in wealth and income inequality.
% average income also increases in equilibrium and now even the
% poorest household have some wealth and minimum income doubled.
Value function iterations: 169 
Wealth distributon mean: 32.5106 
Wealth distributon standard deviation: 18.9376 
Wealth distributon median: 31.9660 
Top 0.1% share of wealth: 0.2923% 
Top 1% share of wealth: 2.6125% 
Top 5% share of wealth: 10.9934% 
Top 10% share of wealth: 20.0534% 
Top 50% share of wealth: 74.5767% 
Maximum wealth: 97.0985 
24
Minimum wealth: 0.1501 
GINI coefficient for wealth: 0.3331 
Income distributon mean: 3.3390 
Income distributon standard deviation: 2.4964 
Income distributon median: 2.8019 
Top 0.1% share of income: 0.4608% 
Top 1% share of income: 4.4347% 
Top 5% share of income: 17.2820% 
Top 10% share of income: 27.5929% 
Top 50% share of income: 73.3018% 
Maximum income: 15.4895 
Minimum income: 0.2143 
GINI coefficient for income: 0.3531 
25
26
27
28
29
Published with MATLAB® R2020b
30