Skip to content

Commit a095d4e

Browse files
robertdkirkbyrobertdkirkby
committed
Residual asset
1 parent f8ac0a9 commit a095d4e

File tree

4 files changed

+986
-0
lines changed

4 files changed

+986
-0
lines changed

‎StationaryDist/FHorz/ResidAsset/StationaryDist_FHorz_Case1_ResidAsset.m

Lines changed: 179 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,179 @@
1+
function StationaryDist=StationaryDist_FHorz_Case1_ResidAsset(jequaloneDist,AgeWeightParamNames,Policy,n_d,n_a,n_z,N_j,pi_z,Parameters,simoptions)
2+
3+
%% Check for the age weights parameter, and make sure it is a row vector
4+
if size(Parameters.(AgeWeightParamNames{1}),2)==1 % Seems like column vector
5+
Parameters.(AgeWeightParamNames{1})=Parameters.(AgeWeightParamNames{1})';
6+
% Note: assumed there is only one AgeWeightParamNames
7+
end
8+
9+
%% Check that the age one distribution is of mass one
10+
if abs(sum(jequaloneDist(:))-1)>10^(-9)
11+
error('The jequaloneDist must be of mass one')
12+
end
13+
14+
%% Setup related to experience asset
15+
if isfield(simoptions,'rprimeFn')
16+
rprimeFn=simoptions.rprimeFn;
17+
else
18+
error('To use an residual asset you must define simoptions.rprimeFn')
19+
end
20+
if ~isfield(simoptions,'a_grid')
21+
error('To use an residual asset you must define simoptions.a_grid')
22+
end
23+
24+
n_r=n_a(end);
25+
n_a=n_a(1:end-1);
26+
27+
a_grid=gpuArray(simoptions.a_grid(1:sum(n_a)));
28+
r_grid=gpuArray(simoptions.a_grid(sum(n_a)+1:end));
29+
30+
% rprimeFnParamNames in same fashion
31+
if n_d(1)==0
32+
l_d=0;
33+
else
34+
l_d=length(n_d);
35+
end
36+
l_a=length(n_a);
37+
l_z=length(n_z);
38+
temp=getAnonymousFnInputNames(rprimeFn);
39+
if length(temp)>(l_d+l_a+l_a+l_z)
40+
rprimeFnParamNames={temp{l_d+l_a+l_a+l_z+1:end}}; % the first inputs will always be (d2,a2)
41+
else
42+
rprimeFnParamNames={};
43+
end
44+
45+
rprimeFnParamNames
46+
47+
%%
48+
if isfield(simoptions,'n_e')
49+
if n_z(1)==0
50+
error('Not yet implemented n_z=0 with n_e and residualasset, email me and I will do it (or you can just pretend by using n_z=1 and pi_z=1, not using the value of z anywhere)')
51+
else
52+
% StationaryDist=StationaryDist_FHorz_Case1_ResidAsset_e(jequaloneDist,AgeWeightParamNames,Policy,n_d1,n_d2,n_a,n_z,N_j,pi_z,rprimeFn,Parameters,simoptions);
53+
end
54+
return
55+
end
56+
57+
if n_z(1)==0
58+
error('Not yet implemented n_z=0 with residualasset, (you can just pretend by using n_z=1 and pi_z=1, not using the value of z anywhere)')
59+
end
60+
61+
N_a=prod(n_a);
62+
N_r=prod(n_r);
63+
N_z=prod(n_z);
64+
65+
if exist('simoptions','var')==0
66+
simoptions.nsims=10^4;
67+
simoptions.parallel=3-(gpuDeviceCount>0); % 3 (sparse) if cpu, 2 if gpu
68+
simoptions.verbose=0;
69+
try
70+
PoolDetails=gcp;
71+
simoptions.ncores=PoolDetails.NumWorkers;
72+
catch
73+
simoptions.ncores=1;
74+
end
75+
simoptions.iterate=1;
76+
simoptions.tolerance=10^(-9);
77+
simoptions.outputkron=0; % If 1 then leave output in Kron form
78+
else
79+
%Check simoptions for missing fields, if there are some fill them with
80+
%the defaults
81+
if isfield(simoptions,'tolerance')==0
82+
simoptions.tolerance=10^(-9);
83+
end
84+
if isfield(simoptions,'nsims')==0
85+
simoptions.nsims=10^4;
86+
end
87+
if isfield(simoptions,'parallel')==0
88+
simoptions.parallel=3-(gpuDeviceCount>0); % 3 (sparse) if cpu, 2 if gpu
89+
end
90+
if isfield(simoptions,'verbose')==0
91+
simoptions.verbose=0;
92+
end
93+
if isfield(simoptions,'ncores')==0
94+
try
95+
PoolDetails=gcp;
96+
simoptions.ncores=PoolDetails.NumWorkers;
97+
catch
98+
simoptions.ncores=1;
99+
end
100+
end
101+
if isfield(simoptions,'iterate')==0
102+
simoptions.iterate=1;
103+
end
104+
if isfield(simoptions,'ExogShockFn') % If using ExogShockFn then figure out the parameter names
105+
simoptions.ExogShockFnParamNames=getAnonymousFnInputNames(simoptions.ExogShockFn);
106+
end
107+
if isfield(simoptions,'outputkron')==0
108+
simoptions.outputkron=0; % If 1 then leave output in Kron form
109+
end
110+
end
111+
112+
jequaloneDistKron=reshape(jequaloneDist,[N_a*N_r*N_z,1]);
113+
if simoptions.parallel~=2 && simoptions.parallel~=4
114+
Policy=gather(Policy);
115+
jequaloneDistKron=gather(jequaloneDistKron);
116+
pi_z=gather(pi_z);
117+
end
118+
119+
% Get policy for aprime, then get policy for rprime, then combine (all just in terms of current state)
120+
Policy=reshape(Policy,[size(Policy,1),N_a*N_r,N_z,N_j]);
121+
122+
if l_a==1
123+
Policy_aprime=shiftdim(Policy(l_d+1,:,:,:),1);
124+
elseif l_a==2
125+
Policy_aprime=shiftdim(Policy(l_d+1,:,:,:)+n_a(1)*(Policy(l_d+2,:,:,:)-1),1);
126+
elseif l_a==3
127+
Policy_aprime=shiftdim(Policy(l_d+1,:,:,:)+n_a(1)*(Policy(l_d+2,:,:,:)-1)+n_a(1)*n_a(2)*(Policy(l_d+3,:,:,:)-1),1);
128+
elseif l_a==4
129+
Policy_aprime=shiftdim(Policy(l_d+1,:,:,:)+n_a(1)*(Policy(l_d+2,:,:,:)-1)+n_a(1)*n_a(2)*(Policy(l_d+3,:,:,:)-1)+n_a(1)*n_a(2)*n_a(3)*(Policy(l_d+4,:,:,:)-1),1);
130+
end
131+
132+
Policy_rprime=zeros(N_a*N_r,N_z,N_j,'gpuArray'); % the lower grid point
133+
PolicyProbs=zeros(N_a*N_r,N_z,N_j,2,'gpuArray'); % The fourth dimension is lower/upper grid point
134+
for jj=1:N_j
135+
rprimeFnParamsVec=CreateVectorFromParams(Parameters, rprimeFnParamNames,jj);
136+
[rprimeIndexes, rprimeProbs]=CreaterprimePolicyResidualAsset_Case1(Policy(:,:,:,jj),rprimeFn, n_d, n_a, n_r, n_z, gpuArray(simoptions.d_grid), a_grid, r_grid, gpuArray(simoptions.z_grid), rprimeFnParamsVec);
137+
% rprimeIndexes is [N_a*N_r,N_z], rprimeProbs is [N_a*N_r,N_z]
138+
Policy_rprime(:,:,jj)=rprimeIndexes;
139+
PolicyProbs(:,:,jj,1)=rprimeProbs;
140+
PolicyProbs(:,:,jj,2)=1-rprimeProbs;
141+
142+
if jj==30 % DEBUG
143+
max(rprimeIndexes)
144+
rprimeIndexes(end-5:end)
145+
rprimeProbs(end-5:end)
146+
end
147+
end
148+
149+
max(max(max(Policy_aprime(:,:,30))))
150+
151+
% size(Policy_aprime)
152+
% size(Policy_rprime)
153+
% max(max(max(abs((Policy_aprime-Policy_rprime).*(PolicyProbs(:,:,:,1)>0)))))
154+
% max(max(max(abs(Policy_aprime-shiftdim(Policy(1,:,:,:),1)))))
155+
156+
Policy_arprime=zeros(N_a*N_r,N_z,N_j,2,'gpuArray');
157+
Policy_arprime(:,:,:,1)=Policy_aprime+N_a*(Policy_rprime-1);
158+
Policy_arprime(:,:,:,2)=Policy_aprime+N_a*(Policy_rprime+1-1);
159+
160+
if simoptions.iterate==0
161+
PolicyProbs=gather(PolicyProbs); % simulation is always with cpu
162+
Policy_arprime=gather(Policy_arprime);
163+
if simoptions.parallel>=3
164+
% Sparse matrix is not relevant for the simulation methods, only for iteration method
165+
simoptions.parallel=2; % will simulate on parallel cpu, then transfer solution to gpu
166+
end
167+
StationaryDistKron=StationaryDist_FHorz_Case1_Simulation_TwoProbs_raw(jequaloneDistKron,AgeWeightParamNames,Policy_arprime,PolicyProbs,N_a*N_r,N_z,N_j,pi_z, Parameters, simoptions);
168+
elseif simoptions.iterate==1
169+
StationaryDistKron=StationaryDist_FHorz_Case1_Iteration_TwoProbs_raw(jequaloneDistKron,AgeWeightParamNames,Policy_arprime,PolicyProbs,N_a*N_r,N_z,N_j,pi_z,Parameters,simoptions); % zero is n_d, because we already converted Policy to only contain aprime
170+
end
171+
172+
if simoptions.outputkron==0
173+
StationaryDist=reshape(StationaryDistKron,[n_a,n_r,n_z,N_j]);
174+
else
175+
% If 1 then leave output in Kron form
176+
StationaryDist=reshape(StationaryDistKron,[N_a,N_r,N_z,N_j]);
177+
end
178+
179+
end

‎ValueFnIter/FHorz/ResidualAsset/ValueFnIter_Case1_FHorz_ResidAsset.m

Lines changed: 105 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,105 @@
1+
function [V,Policy]=ValueFnIter_Case1_FHorz_ResidAsset(n_d,n_a1,n_r,n_z, N_j, d_grid, a1_grid, r_grid, z_grid, pi_z, ReturnFn, Parameters, DiscountFactorParamNames, ReturnFnParamNames, vfoptions)
2+
% vfoptions are already set by ValueFnIter_Case1_FHorz()
3+
if vfoptions.parallel~=2
4+
error('Can only use experience asset with parallel=2 (gpu)')
5+
end
6+
7+
if isfield(vfoptions,'rprimeFn')
8+
rprimeFn=vfoptions.rprimeFn;
9+
else
10+
error('To use an residual asset you must define vfoptions.rprimeFn')
11+
end
12+
13+
if prod(n_d)==0
14+
l_d=0;
15+
else
16+
l_d=length(n_d);
17+
end
18+
l_a1=length(n_a1);
19+
if prod(n_z)==0
20+
l_z=0;
21+
else
22+
l_z=length(n_z);
23+
end
24+
% rprimeFnParamNames in same fashion
25+
temp=getAnonymousFnInputNames(rprimeFn);
26+
if length(temp)>(l_d+l_a1+l_a1+l_z)
27+
rprimeFnParamNames={temp{l_d+l_a1+l_a1+l_z+1:end}}; % the first inputs will always be (d,a1prime,a1,z)
28+
else
29+
rprimeFnParamNames={};
30+
end
31+
32+
N_z=prod(n_z);
33+
34+
if isfield(vfoptions,'n_e')
35+
if isfield(vfoptions,'e_grid_J')
36+
e_grid=vfoptions.e_grid_J(:,1); % Just a placeholder
37+
else
38+
e_grid=vfoptions.e_grid;
39+
end
40+
if isfield(vfoptions,'pi_e_J')
41+
pi_e=vfoptions.pi_e_J(:,1); % Just a placeholder
42+
else
43+
pi_e=vfoptions.pi_e;
44+
end
45+
if n_d==0
46+
if N_z==0
47+
error('Have not implemented residual assets without at least one exogenous variable [you could fake it adding a single-valued z with pi_z=1]')
48+
else
49+
% [VKron, PolicyKron]=ValueFnIter_Case1_FHorz_ResidAsset_nod_e_raw(n_a1,n_r,n_z, vfoptions.n_e, N_j, a1_grid, r_grid, z_grid, e_grid, pi_z, pi_e, ReturnFn, rprimeFn, Parameters, DiscountFactorParamNames, ReturnFnParamNames, rprimeFnParamNames, vfoptions);
50+
end
51+
else
52+
if N_z==0
53+
error('Have not implemented residual assets without at least one exogenous variable [you could fake it adding a single-valued z with pi_z=1]')
54+
else
55+
% [VKron, PolicyKron]=ValueFnIter_Case1_FHorz_ResidAsset_e_raw(n_d,n_a1,n_r,n_z, vfoptions.n_e, N_j, d_grid, a1_grid, r_grid, z_grid, e_grid, pi_z, pi_e, ReturnFn, rprimeFn, Parameters, DiscountFactorParamNames, ReturnFnParamNames, rprimeFnParamNames, vfoptions);
56+
end
57+
end
58+
else
59+
if n_d==0
60+
if N_z==0
61+
error('Have not implemented residual assets without at least one exogenous variable [you could fake it adding a single-valued z with pi_z=1]')
62+
else
63+
[VKron, PolicyKron]=ValueFnIter_Case1_FHorz_ResidAsset_nod_raw(n_a1,n_r,n_z, N_j, a1_grid, r_grid, z_grid, pi_z, ReturnFn, rprimeFn, Parameters, DiscountFactorParamNames, ReturnFnParamNames, rprimeFnParamNames, vfoptions);
64+
end
65+
else
66+
if N_z==0
67+
error('Have not implemented residual assets without at least one exogenous variable [you could fake it adding a single-valued z with pi_z=1]')
68+
else
69+
[VKron, PolicyKron]=ValueFnIter_Case1_FHorz_ResidAsset_raw(n_d,n_a1,n_r,n_z, N_j, d_grid, a1_grid, r_grid, z_grid, pi_z, ReturnFn, rprimeFn, Parameters, DiscountFactorParamNames, ReturnFnParamNames, rprimeFnParamNames, vfoptions);
70+
end
71+
end
72+
end
73+
74+
75+
%%
76+
if vfoptions.outputkron==0
77+
if n_d(1)==0
78+
n_d=n_a1;
79+
else
80+
n_d=[n_d,n_a1];
81+
end
82+
n_a=[n_a1,n_r];
83+
PolicyKron=reshape(PolicyKron,[prod(n_a),prod(n_z),N_j]);
84+
%Transforming Value Fn and Optimal Policy Indexes matrices back out of Kronecker Form
85+
if isfield(vfoptions,'n_e')
86+
if N_z==0
87+
V=reshape(VKron,[n_a,vfoptions.n_e,N_j]);
88+
Policy=UnKronPolicyIndexes_Case2_FHorz(PolicyKron, n_d, n_a, vfoptions.n_e, N_j, vfoptions); % Treat e as z (because no z)
89+
else
90+
V=reshape(VKron,[n_a,n_z,vfoptions.n_e,N_j]);
91+
Policy=UnKronPolicyIndexes_Case2_FHorz_e(PolicyKron, n_d, n_a, n_z, vfoptions.n_e, N_j, vfoptions);
92+
end
93+
else
94+
V=reshape(VKron,[n_a,n_z,N_j]);
95+
Policy=UnKronPolicyIndexes_Case2_FHorz(PolicyKron, n_d, n_a, n_z, N_j, vfoptions);
96+
end
97+
else
98+
V=VKron;
99+
Policy=PolicyKron;
100+
end
101+
102+
103+
end
104+
105+

0 commit comments

Comments
 (0)