%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % Neoclassical growth model using value function iteration % % plots non-linear and linearized policy for % % CRRA or CARA and plots time series % % % % Florian Scheuer, 2006 % % note: uses f.m and u.m functions % %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % housekeeping clear global alpha delta beta sigma CARA %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % parameter values % suggestion: play with the parameters (especially sigma, alpha and delta) % and note how it affects the speed of convergence to the steady state %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% alpha = .33; delta = 1; rho =.03; beta = 1/(1+rho); sigma = 4; % if you want to use the CARA specification set sigma='exp' CARA = 1; % set capital state space interval as fractions of steady state klow = .5; khigh = 1; % grid size and vector of ones gridn = 1000; e = ones(gridn,1); % compute steady state values kss = ((1/beta - 1 + delta)/alpha)^(-1/(1 - alpha)); css = f(kss) - kss; % kgr = ((delta)/alpha)^(-1/(1 - alpha)); % compute golden rule capital % create grid for capital klow = kss*klow; khigh = kss*khigh; k = (klow: (khigh-klow)/(gridn-2) : khigh)'; k = sort([k ; kss]); % add ss capital level % consumption matrix (today's k in column, tomorrow's in rows) C = e*f(k') - k*e' ; % make negative consumption unattractive if sigma ~= 'exp' C = max(C , 0); % take out negative consumption (or else u(c) would not be defined) U = u(C) + ~C*(-1e50); % add penalty for zero consumption else U = u(C); end % initialize iteration crit = 1; iter = 1; v = u(f(k) - k)/(1-beta); % initial guess (value of keeping k and c constant) %v =0*v; % alternative initial guess %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % iterate on Bellman equation % %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% while crit > 1e-20 vnew = max( U + beta*v(:,iter)*e' )'; % compute normalized criterion (scale is in consumption per period) if sigma ~= 'exp' critnew = (1-beta)*max(abs(vnew - v(:,iter)))/css^(-sigma+1); else critnew = (1-beta)*max(abs(vnew - v(:,iter)))/css/exp(-CARA*css); end % display iteration information % note: contraction mapping implies that critnew/crit < beta [iter , critnew, critnew/crit ] crit = critnew; iter = iter +1; v(: , iter) = vnew; end % compute policy function [vnew, policykindex]= max( U + beta*v(:,iter)*e'); policyk = k(policykindex); figure(1); plot(k,[k policyk]); grid %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % plot time series of capital evolution starting % at k0 half the value of the steady state ktime(1) = min(find(k > kss/2)); T = 50; for t=1:T-1 ; ktime(t+1) = policykindex(ktime(t)); end ktime = k(ktime); figure(2); plot( 1:T, ktime/kss); grid title('evolution of capital (as fraction of steady state') %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % solves linearized policy rule (based on linearized EE) ddf = alpha*(alpha-1)*kss^(alpha-2);df = 1/beta; ddUdU = -sigma/css; if sigma == 'exp' ; ddUdU = -CARA; end a = 1; b = -[ 1 + 1/beta + ddf/df/ddUdU ]; c = 1/beta; % compute roots lamda1 = [- b + sqrt(b^2 - 4*a*c)]/a/2; lamda2 = [- b - sqrt(b^2 - 4*a*c)]/a/2; lamda = [lamda1, lamda2]; [bla index] = min(abs(lamda)); % check root with lowest ABSOLUTE value lamda = lamda(index); % take that root % compute linear policy rule policyklinear = (k-kss)*lamda + kss; % plot linear policy rule on top of actual policy function figure(1); hold on; plot(k,policyklinear,'r'); grid on; hold off % plot time series for linear policy on top of previous time series ktimelinear(1)=.5*kss; for t=1:T-1 ; ktimelinear(t+1) = lamda*(ktimelinear(t)-kss)+kss; end figure(2);hold on; plot( 1:T, ktimelinear/kss,'r'); hold off