%SIMULATES AN AIYAGARI ECONOMY WITH IDIOSYNCRATIC BROWNIAN MOTION %GALO NUNO, based on codes from BENJAMIN MOLL %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; 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 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 K = 3.8; % initial aggregate capital. It is important to guess a value close to the solution for the algorithm to converge 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 in 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; %matrices Bj = zeros(I,I); %auxiliary matrix used to construc matrices A and B AA = zeros(I*J); A = zeros(I*J); B = zeros(I*J); 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 Vaf = zeros(I,J); Vab = zeros(I,J); Vzf = zeros(I,J); Vzb = zeros(I,J); Vzz = zeros(I,J); c = zeros(I,J); %CONSTRUCT MATRIX Aswitch SUMMARIZING EVOLUTION OF z yy = - s2/dz2 - mu/dz; chi = s2/(2*dz2); zeta = mu/dz + s2/(2*dz2); Aswitch = zeros(I*J,I*J); %CASE 2<=j<=J-1 for j=2:J-1 for i=1:I row = (j-1)*I+i; i_j = (j-1)*I+i; %index for (i,j) entry i_jp1 = j*I+i; %index for (i,j+1) entry i_jm1 = (j-2)*I+i; %index for (i,j-1) entry Aswitch(row,i_j) = yy(j); Aswitch(row,i_jm1) = chi(j); Aswitch(row,i_jp1) = zeta(j); end end %CASE j=1 j=1; for i=1:I row = (j-1)*I+i; i_j = (j-1)*I+i; %index for (i,j) entry i_jp1 = j*I+i; %index for (i,j+1) entry Aswitch(row,i_j) = yy(j) + chi(j); Aswitch(row,i_jp1) = zeta(j); end %CASE j=J j=J; for i=1:I row = (j-1)*I+i; i_j = (j-1)*I+i; %index for (i,j) entry i_jp1 = j*I+i; %index for (i,j+1) entry i_jm1 = (j-2)*I+i; %index for (i,j-1) entry Aswitch(row,i_j) = yy(j) + zeta(j); Aswitch(row,i_jm1) = chi(j); end %---------------------------------------------------- %INITIAL GUESS r = alpha * K^(alpha-1) -delta; %interest rates w = (1-alpha) * K^(alpha); %wages v0 = (w*zz + r.*aa).^(1-ga)/(1-ga)/rho; v = v0; dist = zeros(1,maxit); %----------------------------------------------------- %MAIN LOOP for iter=1:maxitK disp('Main loop iteration') disp(iter) % 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<=amax just in case % backward difference Vab(2:I,:) = (V(2:I,:)-V(1:I-1,:))/da; Vab(1,:) = (w*z + r.*amin).^(-ga); %state constraint boundary condition I_concave = Vab > 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 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); % 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 < 0; %negative drift --> 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); %CONSTRUCT MATRIX A X = - min(sb,0)/da; Y = - max(sf,0)/da + min(sb,0)/da; Z = max(sf,0)/da; for j=1:J for i=2:I-1 AAj(i,i-1) = X(i,j); AAj(i,i) = Y(i,j); AAj(i,i+1) = Z(i,j); end AAj(1,1)=Y(1,j); AAj(1,2) = Z(1,j); AAj(I,I)=Y(I,j); AAj(I,I-1) = X(I,j); %create block-diagonal matrix AA((j-1)*I+1:j*I,(j-1)*I+1:j*I) = AAj; end A = AA + Aswitch; if max(abs(sum(A,2)))>10^(-12) disp('Improper Transition Matrix') break end B = (1/Delta + rho)*eye(I*J) - A; u_stacked = reshape(u,I*J,1); V_stacked = reshape(V,I*J,1); b = u_stacked + V_stacked/Delta; 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)