Skip to content

Commit 1305b6f

Browse files
robertdkirkbyrobertdkirkby
committed
clean
1 parent a28d4f9 commit 1305b6f

File tree

7 files changed

+144
-19
lines changed

7 files changed

+144
-19
lines changed

‎ValueFnIter/InfHorz/DivideConquer/ValueFnIter_DC1_nod_raw.m

Lines changed: 62 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,62 @@
1+
function [VKron, Policy]=ValueFnIter_Case1_NoD_Par2_raw(VKron, n_a, n_z, pi_z, DiscountFactorParamsVec, ReturnMatrix, Howards,Howards2,Tolerance,maxiter) % Verbose, a_grid, z_grid,
2+
%Does pretty much exactly the same as ValueFnIter_Case1, only without any decision variable (n_d=0)
3+
4+
N_a=prod(n_a);
5+
N_z=prod(n_z);
6+
7+
% PolicyIndexes=zeros(N_a,N_z,'gpuArray');
8+
% Ftemp=zeros(N_a,N_z,'gpuArray');
9+
10+
bbb=reshape(shiftdim(pi_z,-1),[1,N_z*N_z]);
11+
ccc=kron(ones(N_a,1,'gpuArray'),bbb);
12+
aaa=reshape(ccc,[N_a*N_z,N_z]);
13+
14+
addindexforaz=gpuArray(N_a*(0:1:N_a-1)'+N_a*N_a*(0:1:N_z-1));
15+
16+
%%
17+
tempcounter=1;
18+
currdist=Inf;
19+
while currdist>Tolerance && tempcounter<=maxiter
20+
VKronold=VKron;
21+
22+
%Calc the condl expectation term (except beta), which depends on z but not on control variables
23+
EV=VKronold.*shiftdim(pi_z',-1); %kron(ones(N_a,1),pi_z(z_c,:));
24+
EV(isnan(EV))=0; %multilications of -Inf with 0 gives NaN, this replaces them with zeros (as the zeros come from the transition probabilites)
25+
EV=sum(EV,2); % sum over z', leaving a singular second dimension
26+
27+
entireRHS=ReturnMatrix+DiscountFactorParamsVec*EV; %aprime by a by z
28+
29+
%Calc the max and it's index
30+
[VKron,PolicyIndexes]=max(entireRHS,[],1);
31+
32+
tempmaxindex=shiftdim(PolicyIndexes,1)+addindexforaz; % aprime index, add the index for a and z
33+
34+
Ftemp=reshape(ReturnMatrix(tempmaxindex),[N_a,N_z]); % keep return function of optimal policy for using in Howards
35+
36+
PolicyIndexes=PolicyIndexes(:); % a by z (this shape is just convenient for Howards)
37+
VKron=shiftdim(VKron,1); % a by z
38+
39+
VKrondist=VKron(:)-VKronold(:);
40+
VKrondist(isnan(VKrondist))=0;
41+
currdist=max(abs(VKrondist));
42+
43+
% Use Howards Policy Fn Iteration Improvement (except for first few and last few iterations, as it is not a good idea there)
44+
if isfinite(currdist) && currdist/Tolerance>10 && tempcounter<Howards2
45+
for Howards_counter=1:Howards
46+
EVKrontemp=VKron(PolicyIndexes,:);
47+
EVKrontemp=EVKrontemp.*aaa;
48+
EVKrontemp(isnan(EVKrontemp))=0;
49+
EVKrontemp=reshape(sum(EVKrontemp,2),[N_a,N_z]);
50+
VKron=Ftemp+DiscountFactorParamsVec*EVKrontemp;
51+
end
52+
end
53+
54+
tempcounter=tempcounter+1;
55+
56+
end
57+
58+
Policy=reshape(PolicyIndexes,[N_a,N_z]);
59+
60+
61+
62+
end

‎ValueFnIter/InfHorz/DivideConquer/ValueFnIter_DC1_raw.m

Lines changed: 67 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,67 @@
1+
function [VKron, Policy]=ValueFnIter_DC1_raw(VKron, n_d,n_a,n_z, pi_z, beta, ReturnMatrix, Howards,Howards2, Tolerance) %Verbose,
2+
3+
N_d=prod(n_d);
4+
N_a=prod(n_a);
5+
N_z=prod(n_z);
6+
7+
% PolicyIndexes=zeros(N_a,N_z,'gpuArray');
8+
% Ftemp=zeros(N_a,N_z,'gpuArray');
9+
10+
bbb=reshape(shiftdim(pi_z,-1),[1,N_z*N_z]);
11+
ccc=kron(ones(N_a,1,'gpuArray'),bbb);
12+
aaa=reshape(ccc,[N_a*N_z,N_z]);
13+
% I suspect but have not yet double-checked that could instead just use
14+
% aaa=kron(ones(N_a,1,'gpuArray'),pi_z);
15+
16+
addindexforaz=gpuArray(shiftdim((0:1:N_a-1)*N_d*N_a+shiftdim((0:1:N_z-1)*N_d*N_a*N_a,-1),1));
17+
18+
%%
19+
tempcounter=1;
20+
currdist=Inf;
21+
while currdist>Tolerance
22+
VKronold=VKron;
23+
24+
%Calc the condl expectation term (except beta), which depends on z but not on control variables
25+
EV=VKronold.*shiftdim(pi_z',-1);
26+
EV(isnan(EV))=0; %multilications of -Inf with 0 gives NaN, this replaces them with zeros (as the zeros come from the transition probabilites)
27+
EV=sum(EV,2); % sum over z', leaving a singular second dimension
28+
29+
entireRHS=ReturnMatrix+beta*repelem(EV,N_d,1); % z dimension of EV will autoexpand
30+
31+
%Calc the max and it's index
32+
[VKron,PolicyIndexes]=max(entireRHS,[],1);
33+
34+
PolicyIndexes=shiftdim(PolicyIndexes,1);
35+
36+
tempmaxindex=PolicyIndexes+addindexforaz; % aprime index, add the index for a and z
37+
Ftemp=ReturnMatrix(tempmaxindex); % keep return function of optimal policy for using in Howards
38+
39+
VKron=shiftdim(VKron,1); % a by z
40+
41+
% Update currdist
42+
VKrondist=VKron(:)-VKronold(:);
43+
VKrondist(isnan(VKrondist))=0;
44+
currdist=max(abs(VKrondist));
45+
46+
if isfinite(currdist) && currdist/Tolerance>10 && tempcounter<Howards2 %Use Howards Policy Fn Iteration Improvement
47+
Policy_aprime=ceil(PolicyIndexes(:)/N_d);
48+
for Howards_counter=1:Howards
49+
EVKrontemp=VKron(Policy_aprime,:);
50+
51+
EVKrontemp=EVKrontemp.*aaa;
52+
EVKrontemp(isnan(EVKrontemp))=0;
53+
EVKrontemp=reshape(sum(EVKrontemp,2),[N_a,N_z]);
54+
VKron=Ftemp+beta*EVKrontemp;
55+
end
56+
end
57+
58+
tempcounter=tempcounter+1;
59+
60+
end
61+
62+
63+
Policy=zeros(2,N_a,N_z,'gpuArray'); %NOTE: this is not actually in Kron form
64+
Policy(1,:,:)=shiftdim(rem(PolicyIndexes-1,N_d)+1,-1);
65+
Policy(2,:,:)=shiftdim(ceil(PolicyIndexes/N_d),-1);
66+
67+
end

‎ValueFnIter/InfHorz/DivideConquer/ValueFnIter_Case1_DC2B_nod_raw.m renamed to ‎ValueFnIter/InfHorz/DivideConquer/ValueFnIter_DC2B_nod_raw.m

Lines changed: 2 additions & 2 deletions
Original file line numberDiff line numberDiff line change
@@ -1,4 +1,4 @@
1-
function [V,Policy]=ValueFnIter_Case1_DC2B_nod_raw(V0, n_a, n_z, a_grid, z_gridvals, pi_z, ReturnFn, DiscountFactorParamsVec, ReturnFnParamsVec, vfoptions)
1+
function [V,Policy]=ValueFnIter_DC2B_nod_raw(V0, n_a, n_z, a_grid, z_gridvals, pi_z, ReturnFn, DiscountFactorParamsVec, ReturnFnParamsVec, vfoptions)
22
% DC2B: two endogenous states, divide-and-conquer only on the first endogenous state
33

44
N_a=prod(n_a);
@@ -194,4 +194,4 @@
194194

195195

196196

197-
end
197+
end

‎ValueFnIter/InfHorz/DivideConquer/ValueFnIter_Case1_DC2B_raw.m renamed to ‎ValueFnIter/InfHorz/DivideConquer/ValueFnIter_DC2B_raw.m

Lines changed: 2 additions & 2 deletions
Original file line numberDiff line numberDiff line change
@@ -1,4 +1,4 @@
1-
function [V,Policy]=ValueFnIter_Case1_DC2B_raw(V0, n_d, n_a, n_z, d_grid, a_grid, z_gridvals, pi_z, ReturnFn, DiscountFactorParamsVec, ReturnFnParamsVec, vfoptions)
1+
function [V,Policy]=ValueFnIter_DC2B_raw(V0, n_d, n_a, n_z, d_grid, a_grid, z_gridvals, pi_z, ReturnFn, DiscountFactorParamsVec, ReturnFnParamsVec, vfoptions)
22
% DC2B: two endogenous states, divide-and-conquer only on the first endogenous state
33

44
N_d=prod(n_d);
@@ -214,4 +214,4 @@
214214

215215

216216

217-
end
217+
end

‎ValueFnIter/InfHorz/DivideConquer/ValueFnIter_Case1_DC2_nod_raw.m renamed to ‎ValueFnIter/InfHorz/DivideConquer/ValueFnIter_DC2_nod_raw.m

Lines changed: 2 additions & 2 deletions
Original file line numberDiff line numberDiff line change
@@ -1,4 +1,4 @@
1-
function [V,Policy]=ValueFnIter_Case1_DC2_nod_raw(V0, n_a, n_z, a_grid, z_gridvals, pi_z, ReturnFn, DiscountFactorParamsVec, ReturnFnParamsVec, vfoptions)
1+
function [V,Policy]=ValueFnIter_DC2_nod_raw(V0, n_a, n_z, a_grid, z_gridvals, pi_z, ReturnFn, DiscountFactorParamsVec, ReturnFnParamsVec, vfoptions)
22

33
% DC2: two endogenous states, divide-and-conquer both endogenous states
44

@@ -242,4 +242,4 @@
242242

243243

244244

245-
end
245+
end

‎ValueFnIter/InfHorz/DivideConquer/ValueFnIter_Case1_DC2_raw.m renamed to ‎ValueFnIter/InfHorz/DivideConquer/ValueFnIter_DC2_raw.m

Lines changed: 2 additions & 2 deletions
Original file line numberDiff line numberDiff line change
@@ -1,4 +1,4 @@
1-
function [V,Policy]=ValueFnIter_Case1_DC2_raw(V0, n_d, n_a, n_z, d_grid, a_grid, z_gridvals, pi_z, ReturnFn, DiscountFactorParamsVec, ReturnFnParamsVec, vfoptions)
1+
function [V,Policy]=ValueFnIter_DC2_raw(V0, n_d, n_a, n_z, d_grid, a_grid, z_gridvals, pi_z, ReturnFn, DiscountFactorParamsVec, ReturnFnParamsVec, vfoptions)
22
% Divide-and-conquer with two endogenous states
33

44
N_d=prod(n_d);
@@ -266,4 +266,4 @@
266266

267267

268268

269-
end
269+
end

‎ValueFnIter/ValueFnIter_Case1.m

Lines changed: 7 additions & 11 deletions
Original file line numberDiff line numberDiff line change
@@ -487,9 +487,7 @@
487487

488488

489489
%% Divide-and-conquer
490-
if vfoptions.divideandconquer==1 && length(n_a)==1
491-
warning('vfoptions.divideandconquer=1 has not yet been implemented for length(n_a)=1 in infinite horizon, so just ignoring it (and doing vfoptions.divideandconquer=0)')
492-
elseif vfoptions.divideandconquer==1
490+
if vfoptions.divideandconquer==1
493491
if ~isfield(vfoptions,'level1n')
494492
if length(n_a)==1
495493
vfoptions.level1n=51;
@@ -503,26 +501,24 @@
503501

504502
if prod(n_d)==0
505503
if length(n_a)==1
506-
error('Not yet implemented DC1 no d')
507-
% Not sure I can be bothered with this. Will allow big grids, but new GPUs will do this anyway, and is going to be way way slower.
504+
[V,Policy]=ValueFnIter_DC1_nod_raw(V0, n_a, n_z, a_grid, z_gridvals, pi_z, ReturnFn, DiscountFactorParamsVec, ReturnFnParamsVec, vfoptions);
508505
elseif length(n_a)==2
509506
if vfoptions.level1n(2)==n_a(2) % Don't bother with divide-and-conquer on the second endogenous state
510507
vfoptions.level1n=vfoptions.level1n(1); % Only first one is relevant for DC2B
511-
[V,Policy]=ValueFnIter_Case1_DC2B_nod_raw(V0, n_a, n_z, a_grid, z_gridvals, pi_z, ReturnFn, DiscountFactorParamsVec, ReturnFnParamsVec, vfoptions);
508+
[V,Policy]=ValueFnIter_DC2B_nod_raw(V0, n_a, n_z, a_grid, z_gridvals, pi_z, ReturnFn, DiscountFactorParamsVec, ReturnFnParamsVec, vfoptions);
512509
else
513-
[V,Policy]=ValueFnIter_Case1_DC2_nod_raw(V0, n_a, n_z, a_grid, z_gridvals, pi_z, ReturnFn, DiscountFactorParamsVec, ReturnFnParamsVec, vfoptions);
510+
[V,Policy]=ValueFnIter_DC2_nod_raw(V0, n_a, n_z, a_grid, z_gridvals, pi_z, ReturnFn, DiscountFactorParamsVec, ReturnFnParamsVec, vfoptions);
514511
end
515512
end
516513
else % N_d
517514
if length(n_a)==1
518-
error('Not yet implemented DC1')
519-
% Not sure I can be bothered with this. Will allow big grids, but new GPUs will do this anyway, and is going to be way way slower.
515+
[V,Policy]=ValueFnIter_DC1_raw(V0, n_d, n_a, n_z, d_grid, a_grid, z_gridvals, pi_z, ReturnFn, DiscountFactorParamsVec, ReturnFnParamsVec, vfoptions);
520516
elseif length(n_a)==2
521517
if vfoptions.level1n(2)==n_a(2) % Don't bother with divide-and-conquer on the second endogenous state
522518
vfoptions.level1n=vfoptions.level1n(1); % Only first one is relevant for DC2B
523-
[V,Policy]=ValueFnIter_Case1_DC2B_raw(V0, n_d, n_a, n_z, d_grid, a_grid, z_gridvals, pi_z, ReturnFn, DiscountFactorParamsVec, ReturnFnParamsVec, vfoptions);
519+
[V,Policy]=ValueFnIter_DC2B_raw(V0, n_d, n_a, n_z, d_grid, a_grid, z_gridvals, pi_z, ReturnFn, DiscountFactorParamsVec, ReturnFnParamsVec, vfoptions);
524520
else
525-
[V,Policy]=ValueFnIter_Case1_DC2_raw(V0, n_d, n_a, n_z, d_grid, a_grid, z_gridvals, pi_z, ReturnFn, DiscountFactorParamsVec, ReturnFnParamsVec, vfoptions);
521+
[V,Policy]=ValueFnIter_DC2_raw(V0, n_d, n_a, n_z, d_grid, a_grid, z_gridvals, pi_z, ReturnFn, DiscountFactorParamsVec, ReturnFnParamsVec, vfoptions);
526522
end
527523
end
528524
end

0 commit comments

Comments
 (0)