-
Notifications
You must be signed in to change notification settings - Fork 12
/
ff_wkz_vf_vecsv.m
360 lines (310 loc) · 15.6 KB
/
ff_wkz_vf_vecsv.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
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
%% Risky + Safe Asset Dyna Prog Two-Step Solution (Optimized-Vectorized)
% *back to <https://fanwangecon.github.io Fan>'s
% <https://fanwangecon.github.io/CodeDynaAsset/ Dynamic Assets Repository>
% Table of Content.*
%%
function result_map = ff_wkz_vf_vecsv(varargin)
%% FF_WKZ_VF_VECSV solve infinite horizon exo shock + endo asset problem
% This program solves the infinite horizon dynamic savings and risky
% capital asset problem with some ar1 shock. This is the efficient vectorized version
% of
% <https://fanwangecon.github.io/CodeDynaAsset/m_akz/solve/html/ff_wkz_vf.html
% ff_wkz_vf>. See that file for more descriptions.
%
% @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
%
% @include
%
% * <https://fanwangecon.github.io/CodeDynaAsset/m_akz/solve/html/ff_wkz_evf.html ff_wkz_evf>
% * <https://fanwangecon.github.io/CodeDynaAsset/m_akz/paramfunc/html/ffs_akz_set_default_param.html ffs_akz_set_default_param>
% * <https://fanwangecon.github.io/CodeDynaAsset/m_akz/paramfunc/html/ffs_akz_get_funcgrid.html ffs_akz_get_funcgrid>
% * <https://fanwangecon.github.io/CodeDynaAsset/m_akz/solvepost/html/ff_akz_vf_post.html ff_akz_vf_post>
%
% @seealso
%
% * concurrent (safe + risky) loop: <https://fanwangecon.github.io/CodeDynaAsset/m_akz/solve/html/ff_akz_vf.html ff_akz_vf>
% * concurrent (safe + risky) vectorized: <https://fanwangecon.github.io/CodeDynaAsset/m_akz/solve/html/ff_akz_vf_vec.html ff_akz_vf_vec>
% * concurrent (safe + risky) optimized-vectorized: <https://fanwangecon.github.io/CodeDynaAsset/m_akz/solve/html/ff_akz_vf_vecsv.html ff_akz_vf_vecsv>
% * two-stage (safe + risky) loop: <https://fanwangecon.github.io/CodeDynaAsset/m_akz/solve/html/ff_wkz_vf.html ff_wkz_vf>
% * two-stage (safe + risky) vectorized: <https://fanwangecon.github.io/CodeDynaAsset/m_akz/solve/html/ff_wkz_vf_vec.html ff_wkz_vf_vec>
% * two-stage (safe + risky) optimized-vectorized: <https://fanwangecon.github.io/CodeDynaAsset/m_akz/solve/html/ff_wkz_vf_vecsv.html ff_wkz_vf_vecsv>
% * two-stage + interpolate (safe + risky) loop: <https://fanwangecon.github.io/CodeDynaAsset/m_akz/solve/html/ff_iwkz_vf.html ff_iwkz_vf>
% * two-stage + interpolate (safe + risky) vectorized: <https://fanwangecon.github.io/CodeDynaAsset/m_akz/solve/html/ff_iwkz_vf_vec.html ff_iwkz_vf_vec>
% * two-stage + interpolate (safe + risky) optimized-vectorized: <https://fanwangecon.github.io/CodeDynaAsset/m_akz/solve/html/ff_iwkz_vf_vecsv.html ff_iwkz_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_akz_set_default_param(it_param_set);
% Note: param_map and support_map can be adjusted here or outside to override defaults
% param_map('it_w_n') = 50;
% param_map('it_z_n') = 15;
% get armt and func map
[armt_map, func_map] = ffs_akz_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_akz_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_wkz_vf_vecsv';
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_w', 'ar_z'});
[ar_w, ar_z] = params_group{:};
params_group = values(armt_map, {'ar_a_meshk', 'ar_k_mesha', 'mt_coh_wkb', 'it_ameshk_n'});
[ar_a_meshk, ar_k_mesha, mt_coh_wkb, it_ameshk_n] = 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, {'fl_r_save', 'fl_r_borr', 'fl_w',...
'it_z_n', 'fl_crra', 'fl_beta', 'fl_c_min'});
[fl_r_save, fl_r_borr, fl_wage, it_z_n, fl_crra, fl_beta, fl_c_min] = 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_meshk),length(ar_z));
mt_val = mt_val_cur - 1;
mt_pol_a = zeros(length(ar_a_meshk),length(ar_z));
mt_pol_a_cur = mt_pol_a - 1;
mt_pol_k = zeros(length(ar_a_meshk),length(ar_z));
mt_pol_k_cur = mt_pol_k - 1;
mt_pol_idx = zeros(length(ar_a_meshk),length(ar_z));
mt_ev_condi_z_max_kp = zeros(length(ar_w),length(ar_z));
mt_ev_condi_z_max_kp_cur = mt_ev_condi_z_max_kp - 1;
% We did not need these in ff_oz_vf or ff_oz_vf_vec
% see
% <https://fanwangecon.github.io/M4Econ/support/speed/partupdate/fs_u_c_partrepeat_main.html
% fs_u_c_partrepeat_main> for why store using cells.
cl_u_c_store = cell([it_z_n, 1]);
cl_c_valid_idx = cell([it_z_n, 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 Second Stage Problem k*(w,z)
% This is the key difference between this function and
% <https://fanwangecon.github.io/CodeDynaAsset/m_akz/paramfunc/html/ffs_akz_set_functions.html
% ffs_akz_set_functions> which solves the two stages jointly
% Across iterations, only some k*(w,z) is changing, we only update
% consumption matrix when k*(w,z) is changing
% note in ff_wkz_evf a = b
bl_input_override = true;
[mt_ev_condi_z_max, ~, mt_ev_condi_z_max_kp, mt_ev_condi_z_max_bp] = ...
ff_wkz_evf(mt_val_cur, param_map, support_map, armt_map, bl_input_override);
%% Find which k choice differ across iterations?
mt_w_kstar_diff_idx = (mt_ev_condi_z_max_kp_cur ~= mt_ev_condi_z_max_kp);
%% Solve First Stage Problem w*(z) given k*(w,z)
% loop 1: over exogenous states
for it_z_i = 1:length(ar_z)
% State Array fixed
ar_coh_z = mt_coh_wkb(:,it_z_i);
% Get 2nd Stage Choice Arrays
% Update rows where opti k given w=k'+b' is changing
ar_w_kstar_diff_idx = mt_w_kstar_diff_idx(:, it_z_i);
ar_w_kstar_z = mt_ev_condi_z_max_kp(ar_w_kstar_diff_idx, it_z_i);
ar_w_astar_z = mt_ev_condi_z_max_bp(ar_w_kstar_diff_idx, it_z_i);
% Consumption Update
% Note that compared to
% <https://fanwangecon.github.io/CodeDynaAsset/m_akz/paramfunc/html/ffs_akz_set_functions.html
% ffs_akz_set_functions> the mt_c here is much smaller the same
% number of columns (states) as in the ffs_akz_set_functions file,
% but the number of rows equal to ar_w length.
mt_c = f_cons(ar_coh_z', ar_w_astar_z, ar_w_kstar_z);
% EVAL current utility: N by N, f_util defined earlier
if (fl_crra == 1)
mt_utility_update = f_util_log(mt_c);
fl_u_neg_c = f_util_log(fl_c_min);
else
mt_utility_update = f_util_crra(mt_c);
fl_u_neg_c = f_util_crra(fl_c_min);
end
% Eliminate Complex Numbers
mt_it_c_valid_idx = (mt_c <= fl_c_min);
mt_utility_update(mt_it_c_valid_idx) = fl_u_neg_c;
% Update Storage
if (it_iter == 1)
cl_u_c_store{it_z_i} = mt_utility_update;
cl_c_valid_idx{it_z_i} = mt_it_c_valid_idx;
else
cl_u_c_store{it_z_i}(ar_w_kstar_diff_idx,:) = mt_utility_update;
cl_c_valid_idx{it_z_i}(ar_w_kstar_diff_idx,:) = mt_it_c_valid_idx;
end
% EVAL add on future utility, N by N + N by 1
ar_evzp_ak_condi_z = mt_ev_condi_z_max(:, it_z_i);
mt_utility = cl_u_c_store{it_z_i} + fl_beta*ar_evzp_ak_condi_z;
% Index update
% using the method below is much faster than index replace
% see <https://fanwangecon.github.io/M4Econ/support/speed/index/fs_subscript.html fs_subscript>
mt_it_c_valid_idx = cl_c_valid_idx{it_z_i};
mt_utility = mt_utility.*(~mt_it_c_valid_idx) + fl_u_neg_c*(mt_it_c_valid_idx);
% Optimization: remember matlab is column major, rows must be
% choices, columns must be states
% <https://en.wikipedia.org/wiki/Row-_and_column-major_order COLUMN-MAJOR>
[ar_opti_val1_z, ar_opti_idx_z] = max(mt_utility);
mt_val(:,it_z_i) = ar_opti_val1_z;
mt_pol_a(:,it_z_i) = mt_ev_condi_z_max_bp(ar_opti_idx_z, it_z_i);
mt_pol_k(:,it_z_i) = mt_ev_condi_z_max_kp(ar_opti_idx_z, it_z_i);
%% Obtain Choice Index for Vectorized/Analytical Distribution Programs
% For deriving distributions using vectorized and semi-analytical
% methods, at convergence, what index do optimal choices correspond
% to in terms of the rows of mt_val and mt_pol_a, and mt_pol_k. For
% the LHS matrixes here, each column a different optimal choice,
% each row a different element of the a_meshk and k_mesha vectors.
% For each column (one optimal choice), which row has that optimal
% choice's k' and b' values. Note the code is shorter than code in
% <https://fanwangecon.github.io/CodeDynaAsset/m_akz/solve/html/ff_wkz_vf_vec.html
% ff_wkz_vf_vec>.
if (it_iter == (it_maxiter_val + 1))
[~, ar_opti_fullakvec_idx_z] = max((ar_a_meshk == mt_pol_a(:,it_z_i)').*(ar_k_mesha == mt_pol_k(:,it_z_i)'));
mt_pol_idx(:,it_z_i) = ar_opti_fullakvec_idx_z;
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) + norm(mt_pol_k - mt_pol_k_cur);
ar_pol_a_perc_change = sum((mt_pol_a ~= mt_pol_a_cur))/(it_ameshk_n);
ar_pol_k_perc_change = sum((mt_pol_k ~= mt_pol_k_cur))/(it_ameshk_n);
mt_pol_perc_change(it_iter, :) = mean([ar_pol_a_perc_change;ar_pol_k_perc_change]);
% Update
mt_val_cur = mt_val;
mt_pol_a_cur = mt_pol_a;
mt_pol_k_cur = mt_pol_k;
mt_ev_condi_z_max_kp_cur = mt_ev_condi_z_max_kp;
% 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); ...
mean(mt_pol_k_cur,1); ...
mt_val_cur(it_ameshk_n,:); ...
mt_pol_a_cur(it_ameshk_n,:); ...
mt_pol_k_cur(it_ameshk_n,:)]);
tb_valpol_iter.Properties.VariableNames = strcat('z', string((1:size(mt_val_cur,2))));
tb_valpol_iter.Properties.RowNames = {'mval', 'map', 'mak', 'Hval', 'Hap', 'Hak'};
disp('mval = mean(mt_val_cur,1), average value over a')
disp('map = mean(mt_pol_a_cur,1), average choice over a')
disp('mkp = mean(mt_pol_k_cur,1), average choice over k')
disp('Hval = mt_val_cur(it_ameshk_n,:), highest a state val')
disp('Hap = mt_pol_a_cur(it_ameshk_n,:), highest a state choice')
disp('mak = mt_pol_k_cur(it_ameshk_n,:), highest k 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
result_map = containers.Map('KeyType','char', 'ValueType','any');
result_map('mt_val') = mt_val;
result_map('mt_pol_idx') = mt_pol_idx;
result_map('cl_mt_pol_coh') = {mt_coh_wkb, zeros(1)};
result_map('cl_mt_pol_a') = {mt_pol_a, zeros(1)};
result_map('cl_mt_pol_k') = {mt_pol_k, zeros(1)};
result_map('cl_mt_pol_c') = {f_cons(mt_coh_wkb, mt_pol_a, mt_pol_k), zeros(1)};
result_map('ar_st_pol_names') = ["cl_mt_pol_coh", "cl_mt_pol_a", "cl_mt_pol_k", "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_akz_vf_post(param_map, support_map, armt_map, func_map, result_map, bl_input_override);
end
end