% this program solves an Income Fluctuation Problem
% requires obj.m function which is the RHS of the Bellman
% equation as a function of c=consumption and x=cash in hand
global sigma beta yl yh R ppv ppf s1 p1
% possible shocks are in s1 and probability in p1
s1 = [.5 1.5];
p1 = [.5;.5 ];
% other specification to experiment with: (taken from L-S book)
%s1 = [0.3012,0.4493,0.6703,1.0000,1.4918,2.2255,3.3201];
%s2 = [0.0410,0.1769,0.4245,0.8757,1.6978,3.1958,5.9253];
%p1 = [0.0063,0.0608,0.2417,0.3823,0.2417,0.0608,0.0063]';
yl=min(s1); yh=max(s1);
% parameters
% make sure beta*(1+r)<1 !!!
sigma=.5; w = 1; ymean=s1*p1;
beta=.95; delta=1/beta-1
r=.02
R=1+r;
xgrid=100;
% grid for x -- exludes zero to avoid utility==-infinity -- lower bound is loweste income realization
% || note: make sure that xhigh is high enough so that assets contained in our grid interval!
% if you iterate and it seems like this is not the case adjust xhigh accordingly
% likewise, don't make it too large because it is wasteful ||
epsilon=1e-3; xhigh= yh*5; xlow = yl;
x=(xlow : (xhigh - xlow)/(xgrid-1) :xhigh)';
% initial "guess"(es)
c = min((r/1+r)*x + yl/(1+r),x);
v_mean = 1/(1-beta)*(r/R*x + ymean).^(1-sigma)/(1-sigma);
v_aut = (r/R*x).^(1-sigma)/(1-sigma) + beta/(1-beta)*(r/R*x*ones(1,length(s1)) + ones(length(x),1)*s1).^(1-sigma)*p1/(1-sigma);
% use this guess...
v = v_mean;
% ... in spline form (for interpolating between grid points)
ppv=spline(x,v);
% iterate until convergence
iter=0; crit=1; c_new =c; v_new=v;
while crit >1e-10
% compute for each point on the grid for x the optimal c
for i=1:xgrid
% perform the minimization -- obj.m is a function for the negative of the RHS of the bellman
% the optimization is performed over c between 0 and x(i); the last argument is to pass the
% current cash available as a parameter, so that obj is looked at as a function of c only by
% the optimizer
if i>1; c_low=c(i-1) ;else; c_low=0;end; %incorporate that we know consumption increases in x
c_low=0;
[c_new(i) bla]=fminbnd('obj', c_low , x(i) , [] ,x(i));
v_new(i) = -bla; % note: obj is the negative of the true objective function
end
% update interation housekeeping (criterion and iteration number)
v_equiv = ((1-sigma)*(1-beta)*v).^(1/(1-sigma));
v_new_equiv = ((1-sigma)*(1-beta)*v_new).^(1/(1-sigma));
[crit bla1]= max(abs( v_equiv - v_new_equiv ));
[crit_c bla2]= max(abs( (c - c_new)./c )); crit_c= 100*(c_new(bla2) - c(bla2))/c(bla2);
crit_percent= [v_new_equiv(bla1) - v_equiv(bla1)]*100/ymean;
iter=iter+1;
disp([crit_percent , crit_c*100, iter]) ; % display iteration statistics in percentage + or -
disp([x(bla1) , x(bla2)]) ; %displays where the action is on x-grid
v=v_new; c=c_new; ppv=spline(x,v); % updates information after iteration
% iterate on c(x) policy if its quite stable (using Howard iteration technique)
% make sure this does speed up convergence and not make convergence problematic
% removing this part if necessary
if abs(crit_c*100) < 1e-4
disp('howard iteration k=100')
for k=1:100;
v=-obj(c,x);
ppv=spline(x,v);
end
end
end
ppc = spline(x,c); % splines policy function for use in simulation.m
% creates a plot of 45 degree line, a' x_high x_low and c
plot(x,[x x-c R*(x-c)+yh R*(x-c)+yl c]); grid