Matlab function [download]
clearvars tic %start timer %Set parameters for Markov approximation of labor shocks rho=.3;%persistence parameter for labor shocks sigma=.2; %standard deviation of labor shocks m=3; %number of standard deviations from mean mul=1; %mean for labor shocks states=7; %number of possible states for labor shocks [Z,Zprob] = tauchen(states,mul,rho,sigma,m); P=Zprob;%transition matrix z=Z';%possible values for labor shocks alpha=.36;%production parameter MU=3;%degree of risk aversion step=.2; %step size between capital grid beta=.96;%discount factor delta=.08;%depreciation lambda=(1-beta)/beta;%time discount rate it=1;%initialize counter warning off all Zinv=Zprob^100; %%%%%%%%%%% Step One %%%%%%%%%% N0=sum(z.*Zinv(1,:));%stationary employment %%%%%%%%%%% Step Two %%%%%%%%%% K0=5.6; %guess at initial capital stock DifK=1; ik=0; while DifK>1.01 %stopping criteria for aggregate capital stock %%%%%%%%%%% Step Three %%%%%%%%%% w1=(1-alpha)*(K0/N0)^alpha;%wage rate based on capital stock r=(alpha)*(K0/N0)^(alpha-1)-delta%interest rate based on capital stock gridi=[-2:step:40]'; % column i asset holdings grid gridj=[-2:step:40]';% column j asset holdings grid v0=zeros(length(gridj),length(z));%initialize value function g0=zeros(length(gridj),length(z));%initialize policy function jgrid=length(gridj); dif=1;%initialize differences dif2=1; ik=ik+1; epsilon=.001; %%%%%%%%%%% Step Four %%%%%%%%%% display('Performing Value Function Iteration') while dif>epsilon*(1-beta);%stopping criteria for value function iteration for j=1:length(z);%loop over possible labor shocks parfor i=1:length(gridi);%loop over current asset states c=(1+r).*gridi(i)+w1*z(j)-gridj;%consumption if c(1)<0; c(1)=inf; c=(c+.1).*(c>=0);%ensuring nonnegative consumption else c=(c+.1).*(c>=0);%ensuring nonnegative consumption end w=(c.^(1-MU)-1)./(1-MU)+beta.*sum(repmat(P(j,:),length(v0),1).*v0,2);%value function wstar(i,j)=max(w);%finding maximum of value function over grid choice for each state jhold(i,j)=find(max(w)==w);%pointer indicating position of maximum aprime gprime(i,j)=gridj(jhold(i,j));%optimal policy end; end dif=max(max(abs(wstar-v0)));%difference between value functions dif1(it)=dif;%storing difference between value functions dif2=max(max(abs(g0-gprime)));%difference between policy functions policydif(it)=dif2;%storing difference between policy functions g0=gprime;%update policy function v0=wstar;%update value function it=it+1;%update iteration counter end its=100;%increase number of iterations over density function iteration each iteration over K %%%%%%%%%%% Step Five and Step Six %%%%%%%%%% display('Iterating on the distribution function') [f,K1,gridj2]=invarfpar(P,gridj,gprime,step,its);%function that calculates invariant distribution Kstore(ik)=K0%storing aggregate capital over capital iterations DifK=max(abs(K1-K0));%calculating difference in aggregate capital phi=.1;%relaxation parameter for updating aggregate capital stock %%%%%%%%%%% Step Seven %%%%%%%%%% K0=phi*K1+(1-phi)*K0;%update aggregate capital stock display('Dif between r and lambda') r-lambda if r-lambda>0 K0=unifrnd(5,9); else K0=K0; end DifKstore(ik)=DifK%storing difference between aggregate capital iterations end K0=K1; w1=(1-alpha)*(K0/N0)^alpha;%wage based on final value of capital stock r=(alpha)*(K0/N0)^(alpha-1)-delta;%interest rate based on final value of capital stock copt=(1+r).*repmat(gridj,1,length(z))+w1.*repmat(z,length(gridj),1)-gprime;%policy function for consumption figure;plot(gridj,copt) title('Policy function for consumption'); xlabel('Current asset level') ylabel('Optimal Consumption'); legend('Low labor shock','High labor shock') figure;plot(gridj,gprime) title('Policy function for next periods asset'); xlabel('Current asset level') ylabel('Optimal next period asset'); legend('Low labor shock','High labor shock') display('Savings Rate Without Full Insurance') 100*delta*alpha/(r+delta) display('Savings Rate With Full Insurance') 100*delta*alpha/(lambda+delta) toc function[fx]=lininterp(fx1,fx2,x1,x2,x) fx=fx1+((fx2-fx1)./(x2-x1)).*(x-x1); end function[f,muk,gridj2]=invarfpar(P,gridj,gprime,step,its); n=size(P,1); k=size(P,1); step1=step/3; gridj2=min(gridj):step1:max(gridj); m=length(gridj2); f1=zeros(m,n); f2=f1; f2(1,1)=1; f0=(1/(m*n))*ones(m,n); diff=1; it=0; nopt=length(gprime); nk=length(gridj2); while it<its it=it+1; for k=1:n for l=1:m k0=gridj2(l); if k0<=min(gridj2) aprime=gprime(1,:); elseif k0>=max(gridj2) aprime=gprime(nopt,:); else ind=sum(gridj<k0); x1=gridj(ind); x2=gridj(ind+1); fx1=gprime(ind,:); fx2=gprime(ind+1,:); aprime=lininterp(fx1,fx2,x1,x2,gridj2(l)); end kmin=sum(aprime<=min(gridj2)); kmax=sum(aprime>=max(gridj2)); if kmin~=0 krep=n-kmin; indc=find((aprime<=min(gridj2))==1); ku=length(indc)+1; f1(1,[indc])=f1(1,[indc])+P(k,[indc]).*f0(l,k); i=sum((repmat(gridj2',1,n)<=repmat(aprime,m,1)))+1; i=i(ku:n); n0=(aprime(ku:n)-gridj2(i-1))./(gridj2(i)-gridj2(i-1)); kuind=ku:n; for jk=1:length(i); f1(i(jk),kuind(jk))=f1(i(jk),kuind(jk))+P(k,kuind(jk)).*n0(jk).*f0(l,k); f1(i(jk)-1,kuind(jk))=f1(i(jk)-1,kuind(jk))+P(k,kuind(jk)).*(1-n0(jk)).*f0(l,k); end elseif kmax~=0 ku=(n-kmax); indc=find((aprime>=max(gridj2))==1); f1(nk,[indc])=f1(nk,[indc])+P(k,[indc]).*f0(l,k); i=sum((repmat(gridj2',1,n)<=repmat(aprime,m,1)))+1; i=i(1:ku); n0=(aprime(1:ku)-gridj2(i-1))./(gridj2(i)-gridj2(i-1)); kuind=1:ku; for jk=1:length(i); f1(i(jk),kuind(jk))=f1(i(jk),kuind(jk))+P(k,kuind(jk)).*n0(jk).*f0(l,k); f1(i(jk)-1,kuind(jk))=f1(i(jk)-1,kuind(jk))+P(k,kuind(jk)).*(1-n0(jk)).*f0(l,k); end else i=sum((repmat(gridj2',1,n)<=repmat(aprime,m,1)))+1; n0=(aprime-gridj2(i-1))./(gridj2(i)-gridj2(i-1)); kuind=1:n; for jk=1:n; f1(i(jk),kuind(jk))=f1(i(jk),kuind(jk))+P(k,kuind(jk)).*n0(jk).*f0(l,k); f1(i(jk)-1,kuind(jk))=f1(i(jk)-1,kuind(jk))+P(k,kuind(jk)).*(1-n0(jk)).*f0(l,k); end end end end f1=f1./sum(sum(f1)); diff=sum(sum(abs(f1-f0))); difstore(it)=diff; f0=f1; diff; end f=f1; muk=sum(sum(repmat(gridj2',1,n).*f1)); end function [Z,Zprob] = tauchen(N,mu,rho,sigma,m) Z = zeros(N,1); Zprob = zeros(N,N); a = (1-rho)*mu; Z(N) = m * sqrt(sigma^2 / (1 - rho^2)); Z(1) = -Z(N); zstep = (Z(N) - Z(1)) / (N - 1); for i=2:(N-1) Z(i) = Z(1) + zstep * (i - 1); end Z = Z + a / (1-rho); for j = 1:N for k = 1:N if k == 1 Zprob(j,k) = cdf_normal((Z(1) - a - rho * Z(j) + zstep / 2) / sigma); elseif k == N Zprob(j,k) = 1 - cdf_normal((Z(N) - a - rho * Z(j) - zstep / 2) / sigma); else Zprob(j,k) = cdf_normal((Z(k) - a - rho * Z(j) + zstep / 2) / sigma) - ... cdf_normal((Z(k) - a - rho * Z(j) - zstep / 2) / sigma); end end end function c = cdf_normal(x) c = 0.5 * erfc(-x/sqrt(2)); end end