2 while loops one after the other
3 visualizaciones (últimos 30 días)
Mostrar comentarios más antiguos
Panagiotis Veneris
el 5 de Mzo. de 2021
Respondida: Shiva Kalyan Diwakaruni
el 8 de Mzo. de 2021
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?
0 comentarios
Respuestas (1)
Shiva Kalyan Diwakaruni
el 8 de Mzo. de 2021
Hi ,
In the first while loop you have used the variable name 'it' for calculating the iterations and in the second while loop you have used 'iter' for calculating the iterations and used variable 'it' in the if mod(it,100) == 0 change the variable name to 'iter' in
if mod(iter,100) == 0 and you can manage to see the iterations (and whether it has converged) for the second loop in the command window
hope it helps.
0 comentarios
Ver también
Categorías
Más información sobre Loops and Conditional Statements en Help Center y File Exchange.
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!