Hi everyone.
I am trying to make some numerical calculations and I need two different ''while'' loops. The first ''while'' loop calculates A,B for exapmle, and the second ''while'' loop calculates C as a function of A,B.
This is the code
P.gamma = 2;
P.beta = 0.98;
Pr.R = 1/P.beta-1e-4;
tol = 1e-15;
maxit = 10000;
P.gridsize = 200;
P.min = 0;
P.max = 20;
P.grid = linspace(P.min.^0.25, P.max.^0.25, P.gridsize).^4;
P.N = 2;
P.e1 = 0.8;
P.e2 = 1.2;
P.endow = [P.e1,P.e2];
pr = zeros(P.N,P.N);
pr(1,1) = 0.5;
pr(2,1) = 0.5;
pr(1,2) = 0.5;
pr(2,2) = 0.5;
% Auxiliary matrices
P.tiledGrid = repmat(P.grid,2,1);
P.tiledEndow = repmat(P.endow',1,P.gridsize);
% Guess for a
G = 10+0.1*P.tiledGrid;
dif = 1;
it=0;
disp('solving for policy rules')
tic
% Start Loop that Iterates on Euler Equation until the fixed point is reached
while dif>tol && it<maxit
it = it+1;
% Create c_t+1
cp = zeros(P.N,P.gridsize);
for i=1:P.N
app = interp1(G(i,:),P.grid,P.tiledGrid(i,:),'linear','extrap');
app(app<0)=0;
cp(i,:) = Pr.R*P.tiledGrid(i,:)+P.tiledEndow(i,:)-app;
end
upcp = cp.^(-P.gamma); % u'(c_t+1)
Eupcp = [];
for i_p = 1:P.N
Eupcp(i_p,:) = pr(:,i_p)' * upcp;
end
upc = P.beta*Pr.R*Eupcp;
c = upc.^(-1/P.gamma);
a = (P.tiledGrid+c-P.tiledEndow)/Pr.R;
% Check for convergence
if mod(it,50) == 0
dif = max(max(abs(a-G)/(abs(a)+abs(G)+tol)));
fprintf('it = %d, dif = %.3f\n',[it,dif]);
end
G = a;
end
toc
ap(1,:) = interp1(G(1,:),P.grid,P.grid,'linear','extrap');
ap(2,:) = interp1(G(2,:),P.grid,P.grid,'linear','extrap');
ap(ap<P.min)=P.min;
c_pr(1,:) = interp1(G(1,:),c(1,:),P.grid,'linear','extrap');
c_pr(2,:) = interp1(G(2,:),c(2,:),P.grid,'linear','extrap');
u = (c_pr.^(1-P.gamma)-1)./(1-P.gamma);
% Guess for V
V = 0*P.tiledGrid;
test = 1;
iter = 0;
disp('solving for value function')
while test>tol && iter<maxit
iter = iter+1;
v_v = zeros(P.N,P.gridsize);
for i = 1:P.N
%Vp(i,:) = interp1(G(i,:),V(i,:),P.tiledGrid(i,:),'linear','extrap');
Vp(i,:) = interp1(P.grid,V(i,:),ap(i,:),'linear','extrap');
end
E_Vp = [];
for i=1:P.N
E_Vp(i,:) = pr(:,i)'*Vp;
end
V_new = u + P.beta*E_Vp;
if mod(it,100) == 0
test = max(max(abs(V-V_new)/(abs(V)+abs(V_new)+tol)));
fprintf('iter = %d, dif = %.12f\n',[iter,test])
end
V = V_new;
end
figure;
subplot(1,2,1);
hold on;
plot(P.grid,c_pr)
xlabel('Current Assets a');
ylabel('Consumption');
legend({'Low endowment','High endowment'});
%legend('High/Low endowment');
grid on;
subplot(1,2,2);
hold on;
plot(P.grid,ap-P.grid) % y-axis a_t+1, x-axis a_t
xlabel('Current Assets a')
ylabel('Savings a΄')
legend('Low endowment','High endowment');
grid on;
figure;
plot(P.grid,V_new);
xlabel('Current Assets a');
ylabel('V_{t}');
legend('Low endowment','High endowment');
title('Value Function');
grid on;
The thing is that while matlab generates a plot for V_new (calculated in the second ''while'' loop), I am not sure whether the second ''while'' loop works properly. And while in the command window i can see after how many iterations the first iterative procedure has converged, i can't see the same for the second ''while'' loop. Any idea whether the code is fine and how i can manage to see the iterations (and whether it has converged) for the second loop in the command window?
