-
Notifications
You must be signed in to change notification settings - Fork 12
/
ff_az_vf.m
303 lines (262 loc) · 11.3 KB
/
ff_az_vf.m
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
%% Solve One Asset Dynamic Programming Problem (Loop)
% *back to <https://fanwangecon.github.io Fan>'s
% <https://fanwangecon.github.io/CodeDynaAsset/ Dynamic Assets Repository>
% Table of Content.*
%%
function result_map = ff_az_vf(varargin)
%% FF_AZ_VF solve infinite horizon exo shock + endo asset problem
% This program solves the infinite horizon dynamic single asset and single
% shock problem with loops. It is useful to have a version of code that is
% looped for easy debugging. This is the standard dynamic exogenous
% incomplete savings problem.
%
% See
% <https://fanwangecon.github.io/CodeDynaAsset/m_abz/solve/html/ff_abz_vf.html
% ff_abz_vf> for the version of the problem that accommodates both borrowing
% and savings.
%
% @param param_map container parameter container
%
% @param support_map container support container
%
% @param armt_map container container with states, choices and shocks
% grids that are inputs for grid based solution algorithm
%
% @param func_map container container with function handles for
% consumption cash-on-hand etc.
%
% @return result_map container contains policy function matrix, value
% function matrix, iteration results, and policy function, value function
% and iteration results tables.
%
% keys included in result_map:
%
% * mt_val matrix states_n by shock_n matrix of converged value function grid
% * mt_pol_a matrix states_n by shock_n matrix of converged policy function grid
% * ar_val_diff_norm array if bl_post = true it_iter_last by 1 val function
% difference between iteration
% * ar_pol_diff_norm array if bl_post = true it_iter_last by 1 policy
% function difference between iterations
% * mt_pol_perc_change matrix if bl_post = true it_iter_last by shock_n the
% proportion of grid points at which policy function changed between
% current and last iteration for each element of shock
%
% @example
%
% % Get Default Parameters
% it_param_set = 2;
% [param_map, support_map] = ffs_abz_set_default_param(it_param_set);
% % Change Keys in param_map
% param_map('it_a_n') = 50;
% param_map('it_z_n') = 5;
% param_map('fl_a_max') = 100;
% param_map('fl_w') = 1.3;
% % Change Keys support_map
% support_map('bl_display') = false;
% support_map('bl_post') = true;
% support_map('bl_display_final') = false;
% % Call Program with external parameters that override defaults
% % Note this program works very slowly if the grid sizes are too large
% ff_az_vf(param_map, support_map);
%
% @include
%
% * <https://fanwangecon.github.io/CodeDynaAsset/m_az/paramfunc/html/ffs_az_set_default_param.html ffs_az_set_default_param>
% * <https://fanwangecon.github.io/CodeDynaAsset/m_az/paramfunc/html/ffs_az_get_funcgrid.html ffs_az_get_funcgrid>
% * <https://fanwangecon.github.io/CodeDynaAsset/m_az/solvepost/html/ff_az_vf_post.html ff_az_vf_post>
%
% @seealso
%
% * save loop: <https://fanwangecon.github.io/CodeDynaAsset/m_az/solve/html/ff_az_vf.html ff_az_vf>
% * save vectorized: <https://fanwangecon.github.io/CodeDynaAsset/m_az/solve/html/ff_az_vf_vec.html ff_az_vf_vec>
% * save optimized-vectorized: <https://fanwangecon.github.io/CodeDynaAsset/m_az/solve/html/ff_az_vf_vecsv.html ff_az_vf_vecsv>
% * save + borr loop: <https://fanwangecon.github.io/CodeDynaAsset/m_abz/solve/html/ff_abz_vf.html ff_abz_vf>
% * save + borr vectorized: <https://fanwangecon.github.io/CodeDynaAsset/m_abz/solve/html/ff_abz_vf_vec.html ff_abz_vf_vec>
% * save + borr optimized-vectorized: <https://fanwangecon.github.io/CodeDynaAsset/m_abz/solve/html/ff_abz_vf_vecsv.html ff_abz_vf_vecsv>
%
%% Default
% * it_param_set = 1: quick test
% * it_param_set = 2: benchmark run
% * it_param_set = 3: benchmark profile
% * it_param_set = 4: press publish button
it_param_set = 4;
bl_input_override = true;
[param_map, support_map] = ffs_az_set_default_param(it_param_set);
% Note: param_map and support_map can be adjusted here or outside to override defaults
param_map('it_a_n') = 75;
param_map('it_z_n') = 15;
[armt_map, func_map] = ffs_az_get_funcgrid(param_map, support_map, bl_input_override); % 1 for override
default_params = {param_map support_map armt_map func_map};
%% Parse Parameters 1
% if varargin only has param_map and support_map,
params_len = length(varargin);
[default_params{1:params_len}] = varargin{:};
param_map = [param_map; default_params{1}];
support_map = [support_map; default_params{2}];
if params_len >= 1 && params_len <= 2
% If override param_map, re-generate armt and func if they are not
% provided
bl_input_override = true;
[armt_map, func_map] = ffs_az_get_funcgrid(param_map, support_map, bl_input_override);
else
% Override all
armt_map = [armt_map; default_params{3}];
func_map = [func_map; default_params{4}];
end
% append function name
st_func_name = 'ff_az_vf';
support_map('st_profile_name_main') = [st_func_name support_map('st_profile_name_main')];
support_map('st_mat_name_main') = [st_func_name support_map('st_mat_name_main')];
support_map('st_img_name_main') = [st_func_name support_map('st_img_name_main')];
%% Parse Parameters 2
% armt_map
params_group = values(armt_map, {'ar_a', 'mt_z_trans', 'ar_z'});
[ar_a, mt_z_trans, ar_z] = params_group{:};
% func_map
params_group = values(func_map, {'f_util_log', 'f_util_crra', 'f_cons', 'f_coh'});
[f_util_log, f_util_crra, f_cons, f_coh] = params_group{:};
% param_map
params_group = values(param_map, {'it_a_n', 'it_z_n', 'fl_crra', 'fl_beta', 'fl_nan_replace'});
[it_a_n, it_z_n, fl_crra, fl_beta, fl_nan_replace] = params_group{:};
params_group = values(param_map, {'it_maxiter_val', 'fl_tol_val', 'fl_tol_pol', 'it_tol_pol_nochange'});
[it_maxiter_val, fl_tol_val, fl_tol_pol, it_tol_pol_nochange] = params_group{:};
% support_map
params_group = values(support_map, {'bl_profile', 'st_profile_path', ...
'st_profile_prefix', 'st_profile_name_main', 'st_profile_suffix',...
'bl_time', 'bl_display', 'it_display_every', 'bl_post'});
[bl_profile, st_profile_path, ...
st_profile_prefix, st_profile_name_main, st_profile_suffix, ...
bl_time, bl_display, it_display_every, bl_post] = params_group{:};
%% Initialize Output Matrixes
mt_val_cur = zeros(length(ar_a),length(ar_z));
mt_val = mt_val_cur - 1;
mt_pol_a = zeros(length(ar_a),length(ar_z));
mt_pol_a_cur = mt_pol_a - 1;
%% Initialize Convergence Conditions
bl_vfi_continue = true;
it_iter = 0;
ar_val_diff_norm = zeros([it_maxiter_val, 1]);
ar_pol_diff_norm = zeros([it_maxiter_val, 1]);
mt_pol_perc_change = zeros([it_maxiter_val, it_z_n]);
%% Iterate Value Function
% Loop solution with 4 nested loops
%
% # loop 1: over exogenous states
% # loop 2: over endogenous states
% # loop 3: over choices
% # loop 4: add future utility, integration--loop over future shocks
%
% Start Profile
if (bl_profile)
close all;
profile off;
profile on;
end
% Start Timer
if (bl_time)
tic;
end
% Value Function Iteration
while bl_vfi_continue
it_iter = it_iter + 1;
%% Solve Optimization Problem Current Iteration
% loop 1: over exogenous states
for it_z_i = 1:length(ar_z)
fl_z = ar_z(it_z_i);
% loop 2: over endogenous states
for it_a_j = 1:length(ar_a)
fl_a = ar_a(it_a_j);
ar_val_cur = zeros(size(ar_a));
% loop 3: over choices
for it_ap_k = 1:length(ar_a)
fl_ap = ar_a(it_ap_k);
fl_c = f_cons(fl_z, fl_a, fl_ap);
% current utility
if (fl_crra == 1)
ar_val_cur(it_ap_k) = f_util_log(fl_c);
else
ar_val_cur(it_ap_k) = f_util_crra(fl_c);
end
% loop 4: add future utility, integration--loop over future shocks
for it_zp_q = 1:length(ar_z)
ar_val_cur(it_ap_k) = ar_val_cur(it_ap_k) + fl_beta*mt_z_trans(it_z_i,it_zp_q)*mt_val_cur(it_ap_k,it_zp_q);
end
% Replace if negative consumption
if fl_c <= 0
ar_val_cur(it_ap_k) = fl_nan_replace;
end
end
% maximization over loop 3 choices for loop 1+2 states
it_max_lin_idx = find(ar_val_cur == max(ar_val_cur));
mt_val(it_a_j,it_z_i) = ar_val_cur(it_max_lin_idx(1));
mt_pol_a(it_a_j,it_z_i) = ar_a(it_max_lin_idx(1));
end
end
%% Check Tolerance and Continuation
% Difference across iterations
ar_val_diff_norm(it_iter) = norm(mt_val - mt_val_cur);
ar_pol_diff_norm(it_iter) = norm(mt_pol_a - mt_pol_a_cur);
mt_pol_perc_change(it_iter, :) = sum((mt_pol_a ~= mt_pol_a_cur))/(it_a_n);
% Update
mt_val_cur = mt_val;
mt_pol_a_cur = mt_pol_a;
% Print Iteration Results
if (bl_display && (rem(it_iter, it_display_every)==0))
fprintf('VAL it_iter:%d, fl_diff:%d, fl_diff_pol:%d\n', ...
it_iter, ar_val_diff_norm(it_iter), ar_pol_diff_norm(it_iter));
tb_valpol_iter = array2table([mean(mt_val_cur,1); mean(mt_pol_a_cur,1); ...
mt_val_cur(it_a_n,:); mt_pol_a_cur(it_a_n,:)]);
tb_valpol_iter.Properties.VariableNames = strcat('z', string((1:size(mt_val_cur,2))));
tb_valpol_iter.Properties.RowNames = {'mval', 'map', 'Hval', 'Hap'};
disp('mval = mean(mt_val_cur,1), average value over a')
disp('map = mean(mt_pol_a_cur,1), average choice over a')
disp('Hval = mt_val_cur(it_a_n,:), highest a state val')
disp('Hap = mt_pol_a_cur(it_a_n,:), highest a state choice')
disp(tb_valpol_iter);
end
% Continuation Conditions:
% 1. if value function convergence criteria reached
% 2. if policy function variation over iterations is less than
% threshold
if (it_iter == (it_maxiter_val + 1))
bl_vfi_continue = false;
elseif ((it_iter == it_maxiter_val) || ...
(ar_val_diff_norm(it_iter) < fl_tol_val) || ...
(sum(ar_pol_diff_norm(max(1, it_iter-it_tol_pol_nochange):it_iter)) < fl_tol_pol))
% Fix to max, run again to save results if needed
it_iter_last = it_iter;
it_iter = it_maxiter_val;
end
end
% End Timer
if (bl_time)
toc;
end
% End Profile
if (bl_profile)
profile off
profile viewer
st_file_name = [st_profile_prefix st_profile_name_main st_profile_suffix];
profsave(profile('info'), strcat(st_profile_path, st_file_name));
end
%% Process Optimal Choices
% for choices outcomes, store as cell with two elements, first element is
% the y(a,z), outcome given states, the second element will be solved found
% in
% <https://fanwangecon.github.io/CodeDynaAsset/m_az/solve/html/ff_ds_vf.html
% ff_ds_vf> and other distributions files. It stores what are the
% probability mass function of y, along with sorted unique values of y.
result_map = containers.Map('KeyType','char', 'ValueType','any');
result_map('mt_val') = mt_val;
result_map('cl_mt_pol_a') = {mt_pol_a, zeros(1)};
result_map('cl_mt_coh') = {f_coh(ar_z, ar_a'), zeros(1)};
result_map('cl_mt_pol_c') = {f_coh(ar_z, ar_a') - mt_pol_a, zeros(1)};
result_map('ar_st_pol_names') = ["mt_val", "cl_mt_pol_a", "cl_mt_coh", "cl_mt_pol_c"];
if (bl_post)
bl_input_override = true;
result_map('ar_val_diff_norm') = ar_val_diff_norm(1:it_iter_last);
result_map('ar_pol_diff_norm') = ar_pol_diff_norm(1:it_iter_last);
result_map('mt_pol_perc_change') = mt_pol_perc_change(1:it_iter_last, :);
result_map = ff_az_vf_post(param_map, support_map, armt_map, func_map, result_map, bl_input_override);
end
end