clear all; clc; tic; s = 2; r = 0.045; rho = 0.05; w = .1; I=500; amin = -0.02; amax = 1; a = linspace(amin,amax,I)'; da = (amax-amin)/(I-1); maxit=10000; crit = 10^(-6); Delta = 20; dVf = zeros(I,1); dVb = zeros(I,1); c = zeros(I,1); %INITIAL GUESS v0 = (w + r.*a).^(1-s)/(1-s)/rho; v = v0; for n=1:maxit V = v; % forward difference dVf(1:I-1) = (V(2:I)-V(1:I-1))/da; dVf(I) = (w + r.*amax).^(-s); %will never be used, but impose state constraint a<=amax just in case % backward difference dVb(2:I) = (V(2:I)-V(1:I-1))/da; dVb(1) = (w + r.*amin).^(-s); %state constraint boundary condition I_concave = dVb > dVf; %indicator whether value function is concave (problems arise if this is not the case) %consumption and savings with forward difference cf = dVf.^(-1/s); ssf = w + r.*a - cf; %consumption and savings with backward difference cb = dVb.^(-1/s); ssb = w + r.*a - cb; %consumption and derivative of value function at steady state c0 = w + r.*a; dV0 = c0.^(-s); % dV_upwind makes a choice of forward or backward differences based on % the sign of the drift If = ssf > 0; %positive drift --> forward difference Ib = ssb < 0; %negative drift --> backward difference I0 = (1-If-Ib); %at steady state %make sure the right approximations are used at the boundaries %STATE CONSTRAINT: USE BOUNDARY CONDITION UNLESS muf > 0 %Ib(I) = 1; If(I) = 0; dV_Upwind = dVf.*If + dVb.*Ib + dV0.*I0; %important to include third term c = dV_Upwind.^(-1/s); u = c.^(1-s)/(1-s); %CONSTRUCT MATRIX x = -min(ssb,0)/da; y = -max(ssf,0)/da + min(ssb,0)/da; z = max(ssf,0)/da; for i=2:I-1 A(i,i-1) = x(i); A(i,i) = y(i); A(i,i+1) = z(i); end A(1,1)=y(1); A(1,2) = z(1); A(I,I)=y(I); A(I,I-1) = x(I); B = (rho + 1/Delta)*eye(I,I) - A; b = u + V/Delta; V = B\b; %SOLVE SYSTEM OF EQUATIONS Vchange = V - v; v = V; dist(n) = max(abs(Vchange)); if dist(n)