%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % % MODEL 4 -- MARK HUGGETT'S MODEL % % George Hall % Brandeis University % %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% ! rm model4.out diary model4.out; disp('MODEL 4'); disp(''); % % set parameter values % sigma = 1.50; % risk aversion beta = 0.98; % subjective discount factor prob = [ .8 .2; .5 .5]; % prob(i,j) = probability (s(t+1)=sj | s(t) = si) theta = 0.05; % non-interest income if unemployed wage = 1.00; % non-interest income if employed Rstart = 1.021; % initial gross interest rate F = -2.0; % borrowing constraint parameter g = 0.60; % relaxation parameter % % form asset grid % maxast = 8; % maximum value of asset grid minast = -5; % minimum value of asset grid incast = 0.5; % size of asset grid increments nasset = round((maxast-minast)/incast+1); % number of grid points % % loop to find R such that sum(lambda*A) = 0 % liter = 1; maxiter = 50; toler = 0.0001; step = 0.05; R = Rstart; flag = 1; disp('ITERATING ON R'); disp(''); disp(' liter R A newstep'); while (flag ~= 0) & (liter <= maxiter); % % tabulate the utility function such that for zero or negative % consumption utility remains a large negative number so that % such values will never be chosen as utility maximizing % util1=-10000*ones(nasset,nasset); % utility when employed util2=-10000*ones(nasset,nasset); % utility when unemployed for i=1:nasset asset=(i-1)*incast + minast; for j=1:nasset assetp = (j-1)*incast + minast; cons = wage + R*asset - assetp; if assetp >= F & cons > 0; util1(j,i)=(cons)^(1-sigma)/(1-sigma); end; end for j=1:nasset assetp = (j-1)*incast + minast; cons = theta*wage + R*asset - assetp; if assetp >= F & cons > 0; util2(j,i)=(cons)^(1-sigma)/(1-sigma); end; end; end; % % initialize some variables % v = zeros(nasset,2); decis = zeros(nasset,2); test1 = 10; test2 = 10; [rs,cs] = size(util1); % % iterate on Bellman's equation and get the decision % rules and the value function at the optimum % while (test1 ~= 0) | (test2 > .1); for i=1:cs; r1(:,i)=util1(:,i)+beta*(prob(1,1)*v(:,1)+ prob(1,2)*v(:,2)); r2(:,i)=util2(:,i)+beta*(prob(2,1)*v(:,1)+ prob(2,2)*v(:,2)); end; [tv1,tdecis1]=max(r1); [tv2,tdecis2]=max(r2); tdecis=[tdecis1' tdecis2']; tv=[tv1' tv2']; test1=max(any(tdecis-decis)); test2=max(max(abs(tv - v))'); v=tv; decis=tdecis; end; decis=(decis-1)*incast + minast; % % form transition matrix % trans is the transition matrix from state at t (row) % to the state at t+1 (column) % g2=sparse(cs,cs); g1=sparse(cs,cs); for i=1:cs g1(i,tdecis1(i))=1; g2(i,tdecis2(i))=1; end trans=[ prob(1,1)*g1 prob(1,2)*g1; prob(2,1)*g2 prob(2,2)*g2]; trans=trans'; probst = (1/(2*nasset))*ones(2*nasset,1); test = 1; while test > 10^(-8); probst1 = trans*probst; test = max(abs(probst1-probst)); probst = probst1; end; % % vectorize the decision rule to be conformable with probst % calculate new aggregate asset meanA % aa=decis(:); meanA=probst'*aa; % % calculate measure over (k,s) pairs % lambda has same dimensions as decis % lambda=zeros(cs,2); lambda(:)=probst; % % calculate stationary distribution of k % probk=sum(lambda'); % stationary distribution of `captal' probk=probk'; % % form metric and update R % if liter == 1; A = meanA; if meanA > 0.0; step = -step; end; end; Aold = A; Anew = meanA; if sign(Aold) ~= sign(Anew) step = -.5*step; end; disp([ liter R meanA step ]); if abs(step) >= toler; R = R + step; else; flag = 0; end; A = Anew; liter = liter+1; end; % % calculate consumption and expected utility % grid = [ (minast:incast:maxast)' ]; congood = wage*(ones(nasset,1)) + R*grid - grid(tdecis(:,1)); conbad = theta*wage*(ones(nasset,1)) + R*grid - grid(tdecis(:,2)); consum = [congood conbad ]; cons2 = [congood.^2 conbad.^2]; meancon = sum(diag(lambda'*consum)); meancon2 = sum(diag(lambda'*cons2)); varcon = ( meancon2 - meancon^2 ); UTILITY = (consum.^(1-sigma))./(1-sigma); UCEU2 = sum(diag(lambda'*UTILITY)); % % print out results % disp('PARAMETER VALUES'); disp(''); disp(' sigma beta F theta'); disp([ sigma beta F theta]); disp(''); disp('EQUILIBRIUM RESULTS '); disp(''); disp(' R A UCEU meancon varcon'); disp([ R meanA UCEU2 meancon varcon]); % % simulate life histories of the agent % disp('SIMULATING LIFE HISTORY'); a = meanA; % initial level of assets a = 0; n = 5; % number of periods to simulate s0 = 1; % initial state hist = zeros(n-1,2); cons = zeros(n-1,1); invest = zeros(n-1,1); labinc = zeros(n-1,1); grid = [ (minast:incast:maxast)' ]; [chain,state] = markov(prob,n,s0); for i = 1:n-1; hist(i,:) = [ a chain(i) ]; I1 = round((a-minast)/incast) ; I2 = round((a-minast)/incast) + 1; if I1 == 0; I1=1; disp('N.B. I1 = 0'); end; if I2 > nasset; I2 = nasset; disp('N.B. I2 > nasset'); end; weight = (grid(I2,1) - a)/incast; aprime = weight*(decis(I1,chain(i))) + (1-weight)*(decis(I2,chain(i))); if chain(i) == 1; labinc(i) = wage; cons(i) = wage + R*a - aprime; elseif chain(i) == 2; labinc(i) = theta*wage; cons(i) = wage*theta + R*a - aprime; else; disp('something is wrong with chain'); chain end; a = aprime; invest(i) = aprime; end; plot((1:n-1)',invest,(1:n-1)',cons); title('MODEL 4: INVESTMENT AND CONSUMPTION'); print history cov(cons,invest) % % calculate income distribution % income = [ (R*grid + wage) (R*grid + wage*theta) ] ; [ pinc,index ] = sort(income(:)); plambda = lambda(:); plot(pinc,plambda(index)); title('MODEL 4: INCOME DISTRIBUTION'); xlabel('INCOME LEVEL'); ylabel('% OF AGENTS'); print dist %[ grid lambda v ] %save mod4 grid lambda probk %end; diary off