% 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