%%% Example : used to show the basic methods of value function iteration. Junfeng Qiu, CEMA, September, 10, 2007
clear all; %% clear all the variables from the last program
%diary('Q1_3out.txt');
%fid = fopen('Q1_3out.txt','w');
close all; %To close all the figures from the last program
clc %% clear the screen
format long % show 16 digits
% format %% if you want to change back to the default format, just type ``format''
%%%
alpha=0.3; beta=0.6; A=20;
k_low=0.05; %% the lowest value of k in the set of the state variable
k_high=10; %% the highest value of k
k_diff=0.05;
k=[k_low:k_diff:k_high]';
[grid column]=size(k); %% back out the grid
f=zeros(grid,1); %initializing the values of output
for i=1:grid
f(i)=A*(k(i)^alpha); %calculating the output level for each capital input level
end
v=zeros(grid,1); %initializing the values of the value function
tol=0.00001; %setting the tolerance level
diff=tol+1; %initializing the convergence value
vf=ones(grid,1)*(-inf); %initializing function which is for iteration
tv=v; %initializing the values of the iterated value function
k_optimal=zeros(grid,1);
max_loop=100;
v_save=zeros(grid,max_loop); %% this is variable used to save the value of the value function in each iteration
for loop=1:max_loop %% set the maximum loop if the iteration does not converge
v=tv; %replacing the value function by the new guess
v_save(1:end, loop)=v; %% save the value function.
for i=1:grid %starting a loop to get the maximum value function for every today's capital level k(i)
for j=1:grid % given today's k(i), starting a loop to get the maximum value function for tomorrow's capital level k(j)
if f(i)-k(j)<0 %%% this is the case where the capital level is infeasible
vf(j)=-inf; %%% set the value function at very low levels so it will not be chosen
else
vf(j)=log(f(i)-k(j))+beta*v(j); %% calculate the value function given today's capital and tommorow's capital
end
end %end of "for j" loop
[vmax,k_op_index]=max(vf); %% pick the maximum value function, and the index of the optimal k(j)
tv(i)=vmax; %% %replacing each element in the iterated value function
k_optimal(i)=k(k_op_index); %% save the optimal policy function: capital for the next period
end %end of "for i" loop
diff=max(abs(tv-v)); %% repalcing the convergence value
if diff<=tol %% if the value has converged, stop the loop
display('Total number of iterations'); display(loop);
break
end
%display('[k f c_optimal v]')
%display([k f c_optimal v])
end
if loop==max_loop
display('WARNING! Convergence not achieved');
end
k_pine=k_optimal; % optimal capital next period
c_optimal=f-k_pine; % optimal consumption
%%%%%% report the final value and policy functions
format %% use the default number of digits
disp(' the final value and policy functions')
disp(' _____________________________________________________')
disp(' k c k'' v'' ')
disp([k c_optimal k_pine tv])
figure
plot(k, v, 'linewidth', 2.5);
AX=legend('v', 2);
LEG = findobj(AX,'type','text');
set(LEG,'FontSize',21)
set(gca,'FontSize',19)
title('value function')
xlabel('k', 'Fontsize', 19)
ylabel('v','Fontsize',19)
xlim([k_low, k_high]);
print('-depsc2',['example_v.eps']) %% you can use this command to save an eps copy of the figure
figure
plot(k, f, k, c_optimal,'--', k, k_pine, 'linewidth', 2.5);
AX=legend('f(k)', 'optimal consumption','k''', 2);
LEG = findobj(AX,'type','text');
set(LEG,'FontSize',21)
set(gca,'FontSize',19)
title('optimal policy')
xlabel('k', 'Fontsize', 19)
xlim([k_low, k_high]);
figure
plot(k, k_pine,k, k, '--', 'linewidth', 2.5);
AX=legend('Policy function k''','45 degree line', 4);
LEG = findobj(AX,'type','text');
set(LEG,'FontSize',21)
set(gca,'FontSize',19)
title('policy function')
xlabel('k', 'Fontsize', 19)
ylabel('k''', 'Fontsize',19)
xlim([k_low, k_high]);
figure
set(gcf, 'Renderer', 'painters')
plot(k, v_save(:,1), k, v_save(:,2), k, v_save(:, 3),k, v_save(:,6), k, v_save(:,11),k, v_save(:,21), 'linewidth', 2.5);
AX=legend('v0', 'v1', 'v2', 'v5', 'v10', 'v20', 4);
LEG = findobj(AX,'type','text');
set(LEG,'FontSize',21)
set(gca,'FontSize',19)
title('the convergence of the value function')
xlabel('k', 'Fontsize', 19)
print('-depsc2',['v_convergence.eps'])
%%%% compute the analytical solution
V_true=zeros(grid, 1);
k_pine_true=zeros(grid, 1);
for ii=1: grid
V_true(ii)=1/(1-beta)*(log(A*(1-alpha*beta))+beta*alpha/(1-alpha*beta)*log(alpha*beta*A))+alpha/(1-alpha*beta)*log(k(ii));
k_pine_true(ii)=A*beta*alpha*k(ii)^alpha;
end
disp(' _____________________________________________________')
disp(' k k''(numerical) k''(true) v(numerical) v(true)')
disp([k f-c_optimal k_pine_true, tv, V_true ])
figure
plot(k, k_pine,k, k_pine_true, 'linewidth', 2.5);
AX=legend('k''', 'k''_{true}', 2);
LEG = findobj(AX,'type','text');
set(LEG,'FontSize',21)
set(gca,'FontSize',19)
xlabel('k', 'Fontsize', 19)
xlim([k_low, k_high]);
figure
plot(k, v, k, V_true, 'linewidth', 2.5);
AX=legend('v','V_{true}', 2);
LEG = findobj(AX,'type','text');
set(LEG,'FontSize',21)
set(gca,'FontSize',19)
title('value function')
xlabel('k', 'Fontsize', 19)
ylabel('v','Fontsize',19)
xlim([k_low, k_high]);
% fclose(fid); %% close the output file
% diary off
