Finding the general equilibrium

1

I was looking for a faster alternative to fminsearch (the default used by the toolkit) and wanted to test fsolve, but I think it is not implemented among the options heteroagentoptions.fminalgo.

fsolve might be faster in some cases because is based on a trust-region algorithm while fminsearch is a derivative-free method.

With fminsearch the objective function to be minimized is the sum of squared deviations (a scalar)

sqrt(sum(heteroagentoptions.multiGEweights.*(GeneralEqmConditionsVec.^2)))

while fsolve solves a system of equations, so the objective should be simply GeneralEqmConditionsVec=0, I think.

By the way, I realized now that there is multiGEweights. Can I give a different weight to different GE conditions? How do I specifify this? Moreover, does it make sense to give different weights to GE conditions (this is something I do in estimation)?

Thanks!

2

I will look at adding fsolve as an option, Iā€™ll get back to you.

For heteroagentoptions.multiGEweights, internally by default this is a vector of ones. You can change it to any other vector.

Given that all GE conditions should be able to jointly evaluate to zero, the weights are irrelevant (e.g. weights.*zeros(5,1)=0, regardless of the weights themselves).

As you note though, when using the GE commands for calibration, it makes sense to use weights. Specifically, put higher weights on the GE conditions (as it is important that these evaluate to zero) and lower weights on your calibration targets. You might also want to weight some calibration targets higher/lower than other calibration targets.

I donā€™t currently have a way to set weights based on the names of the GE conditions, maybe I should. If you look at the GeneralEqmEqns structure (or whatever name exactly you have given it) you will see the order of the field names, and you should set the weights in this same order.

There is already one way to do GE faster, namely using a shooting-algorithm, see the Appendix ā€œHow to solve faster, or solve more robustlyā€ of Intro to OLG models.

1 Like
3

Thanks for your detailed answer! I will look into the shooting option

4

fsolve can now be accessed setting heteroagentoptions.fminalgo=7

1 Like
5

Great thanks! Iā€™ll test this option in my BueraShin repo

6

I ran the GE with heteroagentoptions.fminalgo=7 but I got this error

Error using fzero>localFirstFcnEvalError (line 733)
FZERO cannot continue because supplied function_handle ==>
@(p)HeteroAgentStationaryEqm_Case1_subfn(p,n_d,n_a,n_z,l_p,pi_z,d_grid,a_grid,z_grid,ReturnFn,FnsToEvaluate,GeneralEqmEqns,Parameters,DiscountFactorParamNames,ReturnFnParamNames,FnsToEvaluateParamNames,GeneralEqmEqnParamNames,GEPriceParamNames,heteroagentoptions,simoptions,vfoptions)
failed with the error below.

Index exceeds the number of array elements. Index must not exceed 1.

Error in fzero (line 231)
localFirstFcnEvalError(FunFcn,FunFcnIn,ME);

Error in HeteroAgentStationaryEqm_Case1 (line 295)
[p_eqm_vec,GeneralEqmConditions]=fzero(GeneralEqmConditionsFnOpt,p0,minoptions);

Error in BueraShin_Fn (line 12)
[p_eqm,~,~]=HeteroAgentStationaryEqm_Case1(n_d, n_a, n_z, 0, pi_z, d_grid, a_grid, z_grid, ReturnFn, FnsToEvaluate, GeneralEqmEqns, Params, DiscountFactorParamNames, , , , GEPriceParamNames,heteroagentoptions, simoptions, vfoptions);

Error in BueraShin2013 (line 150)
[Outputs,GE_cond,Policy_init,StationaryDist_init,AggVars_init,ValuesOnGrid] = BueraShin_Fn(do_GE,Params,n_d,n_a,n_z,pi_z,d_grid,a_grid,z_grid,ReturnFn,FnsToEvaluate,GeneralEqmEqns,DiscountFactorParamNames,GEPriceParamNames,heteroagentoptions,simoptions,vfoptions);

7

Haha, that was stupid! I added an fminalgo=7 to do fsolve() and then just had it do fzero() anyway! Fixed now.

1 Like
8

This time I got this warning:

Current GE prices: 
	r:   0.0476 
	w:   0.1720 
Current aggregate variables: 
	A:   0.7441 
	K:   0.7453 
	L:   0.9466 
	entrepreneur:   0.0533 
	Y:   0.3076 
	extfin:   0.5276 
	earnings:   0.2629 
Current GeneralEqmEqns: 
	capitalmarket:   0.0012 
	labormarket:  -0.0002 
Warning: Trust-region-dogleg algorithm of FSOLVE cannot handle non-square systems; using
Levenberg-Marquardt algorithm instead. 
> In fsolve (line 348)
In HeteroAgentStationaryEqm_Case1 (line 295)
In BueraShin_Fn (line 12)
In BueraShin2013 (line 151)
9

I looked under the hood in HeteroAgentStationaryEqm_Case1_subfn and the reason for the warning is the following:

I have two GE variables (r and w) and two GE conditions, but inside HeteroAgentStationaryEqm_Case1_subfn the code takes the sum of squares of the GE conditions and returns GeneralEqmConditions as a scalar. This is needed for fminsearch (which is a minimizer) but not for fsolve which is a root-finding solver. Now Matlab thinks that the system of equations has two variables but only on equation.

It would be better to have a flag in the code that if fminalgo==7, then the output of HeteroAgentStationaryEqm_Case1_subfn is the vector of the GE conditions. In this way, Matlab uses the Trust-region-dogleg algorithm which is good (in my view)

1 Like
10

Easy fix, there is heteroagentoptions.ouputGEform (for internal purposes only, not something the user should be setting). =0 means ā€œsubfnā€ returns scalar, =1 returns vector. I have updated the GE codes so it =1 when you use fminalgo=7. Have just pushed this :grin:

[As you can probably tell, Iā€™ve never used fsolve(). The shooting algorithms, fminalgo=5, also require returning a vector, hence why the functionality is already there.]

1 Like
11

There is another error

Error using fsolve (line 319)
FSOLVE requires all values returned by functions to be of data type double.

Error in HeteroAgentStationaryEqm_Case1 (line 297)
        [p_eqm_vec,GeneralEqmConditions]=fsolve(GeneralEqmConditionsFnOpt,p0,minoptions);

Error in BueraShin_Fn (line 12)
    [p_eqm,~,~]=HeteroAgentStationaryEqm_Case1(n_d, n_a, n_z, 0, pi_z, d_grid, a_grid, z_grid, ReturnFn, FnsToEvaluate, GeneralEqmEqns, Params, DiscountFactorParamNames, [], [], [], GEPriceParamNames,heteroagentoptions, simoptions, vfoptions);

Error in BueraShin2013 (line 151)
[Outputs,GE_cond,Policy_init,StationaryDist_init,AggVars_init,ValuesOnGrid] = BueraShin_Fn(do_GE,Params,n_d,n_a,n_z,pi_z,d_grid,a_grid,z_grid,ReturnFn,FnsToEvaluate,GeneralEqmEqns,DiscountFactorParamNames,GEPriceParamNames,heteroagentoptions,simoptions,vfoptions);

GeneralEqmConditions is a gpuArray but fsolve accepts only arrays on the GPU (I think). With fminsearch this problem does not occur because GeneralEqmConditions is always a scalar.

Sorry to give you so much trouble :slight_smile:

By the way, I got the idea of using fsolve to solve for the GE after reading a paper by BKM " Exploiting MIT Shocks in Heterogeneous-Agent Economies: The Impulse Response as a Numerical Derivative"

image
image927Ɨ263 40.5 KB

12

Iā€™m away for next few days, will try a ā€œgather()ā€ when I get back, if it is just about returning a value on CPU rather than GPU that will fix it. Will post once I do.

Off-topic, but BKM can easily be implemented building on using VFI Toolkit to solve the (certainty equivalent) transition path [a.k.a. what they refer to as the impulse reponse]. Model simulation is then just a weighted sum of this. I decided not to add it as a feature as it only simulates the aggregate economy, and not the cross section so would fit awkwardly with the rest of the toolkit. [Obviously it is not about to be as fast as taking BKM idea of ā€˜linearization to make everything just weighted IRFsā€™ to the extreme to get SSJ.]

Interesting to see BKM use fsolve().

1 Like
13

Should be using gather() by default when heteroagentoptions.fminalgo=7 [fsovle()]

Hope that works.

1 Like
14

Thank you so much! I will check and write again if there are problems

15

Now it works and I can confirm that fsolve is much faster than fminsearch. This is true for this particular example of course, but I suspect it might hold in general. The reason is that the trust-region method (fsolve) is typically much faster than Nelder-Mead (fminsearch). Nelder-Mead is however more robust since it is a derivative-free method.

In my projects I tend to use fsolve more to find GE, while I use fminsearch to do calibration/estimation. They are both local algorithms though. So for a difficult calibration/estimation, better to use other methods (e.g. simulated annealing or CMA-ES).

1 Like
16

Hmm, maybe I need to run a whole series of tests and if fsolve() is reliably faster for general eqm then I should change the default. Will think about this later in the year if I can find some spare time.

1 Like
17

I have another small request :slight_smile:

Does the VFI toolkit allow for fzero? fzero is very fast since it is based on Brentā€™s rootfinding algorithm and is my preferred choice if I have only one GE variable (fzero works only for one-variable rootfinding problems).

Note: fzero accepts both a scalar or an interval as initial condition. I think the toolkit should do something like

GEPriceParamNames={'r'};
par.r = 0.0317;

and inside HeteroAgentStationaryEqm_Case1

p0(1)=Parameters.(GEPriceParamNames{1});
if numel(GEPriceParamNames)>1
    error('Cannot use fzero if more than one GE variable')
else
 [p_eqm_vec,GeneralEqmConditions]=fzero(GeneralEqmConditionsFnOpt,p0,options_fzero);
end

Maybe the toolkit cannot dealt with the initial condition for r being an interval? But at least could do the scalar case.

18

I tested fzero vs fsolve vs fminsearch on an Aiyagari model and I can confirm that fzero (with scalar initial condition) is much faster (3 or 4 times faster).

So whenever there is only one GE condition and one GE variable, fzero is the best option.

19

Changed all the stationary general eqm commands in toolkit to use the default of fzero when only one general eqm conditions, fminsearch when more than one:

if length(fieldnames(GeneralEqmEqns))==1 
  heteroagentoptions.fminalgo=0; % use fzero
else
  heteroagentoptions.fminalgo=1; % use fminsearch
end

I havenā€™t tried the setting an interval for the general eqm parameters. My feeling is that given models with just one general eqm equation solve fast anyway, there would be very few users who want to try and find ways that make them faster (ways that can be hard-coded, like using fzero instead of fminsearch are good, but ways that require user to actively do things need a higher bar).

Do you feel like the interval gives a meaningful speed gain? Or is it fairly minor?

PS. fminalgo=0 for fzero() already existed, I have just changed the default settings.

1 Like
20

Hi Robert,

Thanks for implementing this change. However, if I run your Aiyagari example, I get an error after a few iterations: fzero sets some strange values for the interest rate which break the code (Technically, we get complex numbers, but this is irrelevant).

The reason is because the function that you feed to fzero is the sum of squared deviation, but it should be simply GeneralEqmConditionsVec, since fzero is a root-finding algorithm, not a minimization algorithm as fminsearch.

So, Iā€™d suggest to modify HeteroAgentStationaryEqm_Case1_subfn to deal with fminalgo=0. I think it should be the same of fminalgo=7 (fsolve, also a rootfinding)

Here is the actual (likely not correct) code:

GeneralEqmConditions=sqrt(sum(heteroagentoptions.multiGEweights.*(GeneralEqmConditionsVec.^2)));

Instead, the code should go here (I think, but please double-check me :smiley:

elseif heteroagentoptions.outputGEform==1 % vector
    GeneralEqmConditions=GeneralEqmConditionsVec;
    if heteroagentoptions.outputgather==1
        GeneralEqmConditions=gather(GeneralEqmConditions);