/
HeteroAgentStationaryEqm_Case2.m
162 lines (145 loc) · 7.55 KB
/
HeteroAgentStationaryEqm_Case2.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
function [p_eqm,p_eqm_index,GeneralEqmConditions]=HeteroAgentStationaryEqm_Case2(n_d, n_a, n_s, n_p, pi_s, d_grid, a_grid, s_grid,Phi_aprimeKron, Case2_Type, ReturnFn, FnsToEvaluate, GeneralEqmEqns, Parameters, DiscountFactorParamNames, ReturnFnParamNames, PhiaprimeParamNames, FnsToEvaluateParamNames, GeneralEqmEqnParamNames, GEPriceParamNames,heteroagentoptions, simoptions, vfoptions)
% If n_p=0 then will use fminsearch to find the general equilibrium (find
% price vector that corresponds to MarketClearance=0). By setting n_p to
% nonzero it is assumend you want to use a grid on prices, which must then
% be passed in heteroagentoptions.p_grid
N_d=prod(n_d);
N_a=prod(n_a);
N_s=prod(n_s);
N_p=prod(n_p);
if isempty(n_p)
N_p=0;
end
l_p=length(GEPriceParamNames);
p_eqm=struct();
%% Check which simoptions and vfoptions have been used, set all others to defaults
if exist('vfoptions','var')==0
%If vfoptions is not given, just use all the defaults
vfoptions.parallel=1+(gpuDeviceCount>0);
%Note that the defaults will be set when we call 'ValueFnIter...'
%commands and the like, so no need to set them here except for a few.
else
%Check vfoptions for missing fields, if there are some fill them with the defaults
if isfield(vfoptions,'parallel')==0
vfoptions.parallel=1+(gpuDeviceCount>0);
end
end
if exist('simoptions','var')==0
simoptions.parallel=1+(gpuDeviceCount>0);
%Note that the defaults will be set when we call 'StationaryDist...'
%commands and the like, so no need to set them here except for a few.
else
%Check simoptions for missing fields, if there are some fill them with the defaults
if isfield(simoptions,'parallel')==0
simoptions.parallel=1+(gpuDeviceCount>0);
end
end
if exist('heteroagentoptions','var')==0
heteroagentoptions.multiGEcriterion=1;
heteroagentoptions.multiGEweights=ones(1,length(GeneralEqmEqns));
heteroagentoptions.toleranceGEprices=10^(-4); % Accuracy of general eqm prices
heteroagentoptions.toleranceGEcondns=10^(-4); % Accuracy of general eqm eqns
heteroagentoptions.verbose=0;
heteroagentoptions.fminalgo=1; % use fminsearch
heteroagentoptions.outputGEform=0; % output of subfn is a scalar
else
if isfield(heteroagentoptions,'multiGEcriterion')==0
heteroagentoptions.multiGEcriterion=1;
end
if isfield(heteroagentoptions,'multiGEweights')==0
heteroagentoptions.multiGEweights=ones(1,length(GeneralEqmEqns));
end
if N_p~=0
if isfield(heteroagentoptions,'pgrid')==0
disp('VFI Toolkit ERROR: you have set n_p to a non-zero value, but not declared heteroagentoptions.pgrid')
end
end
if isfield(heteroagentoptions,'toleranceGEprices')==0
heteroagentoptions.toleranceGEprices=10^(-4); % Accuracy of general eqm prices
end
if isfield(heteroagentoptions,'toleranceGEcondns')==0
heteroagentoptions.toleranceGEcondns=10^(-4); % Accuracy of general eqm prices
end
if isfield(heteroagentoptions,'verbose')==0
heteroagentoptions.verbose=0;
end
if isfield(heteroagentoptions,'fminalgo')==0
heteroagentoptions.fminalgo=1; % use fminsearch
end
if isfield(heteroagentoptions,'outputGEform')==0
heteroagentoptions.outputGEform=0; % output of subfn is a scalar
end
end
%%
zerosinphi_aprimekron=sum(sum(sum(sum(Phi_aprimeKron==0))));
fprintf('If this number is not zero there is an problem with Phi_aprimeKron: %.0f', zerosinphi_aprimekron)
%%
% Check if there is an initial guess for V0
if isfield(vfoptions,'V0')
vfoptions.V0=reshape(vfoptions.V0,[N_a,N_s]);
else
if vfoptions.parallel==2
vfoptions.V0=zeros([N_a,N_s], 'gpuArray');
else
vfoptions.V0=zeros([N_a,N_s]);
end
end
if N_p~=0
[p_eqm_vec,p_eqm_index,GeneralEqmConditions]=HeteroAgentStationaryEqm_Case2_pgrid(n_d, n_a, n_s, n_p, pi_s, d_grid, a_grid, s_grid, Phi_aprimeKron, Case2_Type, ReturnFn, FnsToEvaluate, GeneralEqmEqns, Parameters, DiscountFactorParamNames, ReturnFnParamNames, PhiaprimeParamNames, FnsToEvaluateParamNames, GeneralEqmEqnParamNames, GEPriceParamNames,heteroagentoptions, simoptions, vfoptions);
for ii=1:length(GEPriceParamNames)
p_eqm.(GEPriceParamNames{ii})=p_eqm_vec(ii);
end
return
end
%% Otherwise, use fminsearch to find the general equilibrium
GeneralEqmConditionsFnOpt=@(p) HeteroAgentStationaryEqm_Case2_subfn(p, n_d, n_a, n_s,l_p, pi_s, d_grid, a_grid, s_grid, Phi_aprimeKron, Case2_Type, ReturnFn, FnsToEvaluate, GeneralEqmEqns, Parameters, DiscountFactorParamNames, ReturnFnParamNames, PhiaprimeParamNames, FnsToEvaluateParamNames, GeneralEqmEqnParamNames, GEPriceParamNames,heteroagentoptions, simoptions, vfoptions)
p0=nan(length(GEPriceParamNames),1);
for ii=1:length(GEPriceParamNames)
p0(ii)=Parameters.(GEPriceParamNames{ii});
end
% Choosing algorithm for the optimization problem
% https://au.mathworks.com/help/optim/ug/choosing-the-algorithm.html#bscj42s
if heteroagentoptions.fminalgo==0 % fzero doesn't appear to be a good choice in practice, at least not with it's default settings.
heteroagentoptions.multiGEcriterion=0;
[p_eqm_vec,GeneralEqmConditions]=fzero(GeneralEqmConditionsFnOpt,p0);
elseif heteroagentoptions.fminalgo==1
[p_eqm_vec,GeneralEqmConditions]=fminsearch(GeneralEqmConditionsFnOpt,p0);
elseif heteroagentoptions.fminalgo==2
% Use the optimization toolbox so as to take advantage of automatic differentiation
z=optimvar('z',length(p0));
optimfun=fcn2optimexpr(GeneralEqmConditionsFnOpt, z);
prob = optimproblem("Objective",optimfun);
z0.z=p0;
[sol,GeneralEqmConditions]=solve(prob,z0);
p_eqm_vec=sol.z;
% Note, doesn't really work as automattic differentiation is only for
% supported functions, and the objective here is not a supported function
elseif heteroagentoptions.fminalgo==3
goal=zeros(length(p0),1);
weight=ones(length(p0),1); % I already implement weights via heteroagentoptions
[p_eqm_vec,GeneralEqmConditionsVec] = fgoalattain(GeneralEqmConditionsFnOpt,p0,goal,weight);
GeneralEqmConditions=sum(abs(GeneralEqmConditionsVec));
elseif heteroagentoptions.fminalgo==4 % CMA-ES algorithm (Covariance-Matrix adaptation - Evolutionary Stategy)
% https://en.wikipedia.org/wiki/CMA-ES
% https://cma-es.github.io/
% Code is cmaes.m from: https://cma-es.github.io/cmaes_sourcecode_page.html#matlab
if ~isfield(heteroagentoptions,'insigma')
% insigma: initial coordinate wise standard deviation(s)
heteroagentoptions.insigma=0.3*abs(p0)+0.1*(p0==0); % Set standard deviation to 30% of the initial parameter value itself (cannot input zero, so add 0.1 to any zeros)
end
if ~isfield(heteroagentoptions,'inopts')
% inopts: options struct, see defopts below
heteroagentoptions.inopts=[];
end
% varargin (unused): arguments passed to objective function
if heteroagentoptions.verbose==1
disp('VFI Toolkit is using the CMA-ES algorithm, consider giving a cite to: Hansen, N. and S. Kern (2004). Evaluating the CMA Evolution Strategy on Multimodal Test Functions' )
end
% This is a minor edit of cmaes, because I want to use 'GeneralEqmConditionsFnOpt' as a function_handle, but the original cmaes code only allows for 'GeneralEqmConditionsFnOpt' as a string
[p_eqm_vec,GeneralEqmConditions,counteval,stopflag,out,bestever] = cmaes_vfitoolkit(GeneralEqmConditionsFnOpt,p0,heteroagentoptions.insigma,heteroagentoptions.inopts); % ,varargin);
end
p_eqm_index=nan; % If not using p_grid then this is irrelevant/useless
for ii=1:length(GEPriceParamNames)
p_eqm.(GEPriceParamNames{ii})=p_eqm_vec(ii);
end
end