General equilibrium conditions, weights

1

Question about GE conditions. Sometimes we would like to give more weight to some conditions relative to others. Obviously, it does not matter theoretically as long as the number of GE parameters is equal to the number of GE conditions. In practice, however, when using a minimization routine like fminsearch, you can make convergence faster by giving a larger weight to some conditions that are more difficult to clear.

We have seen that in HeteroAgentStationaryEqm_Case1_subfn.m there is a variable called heteroagentoptions.multiGEweights. How does it work? Can the user set it?

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

In our code we have two GE conditions:

GeneralEqmEqns.CapitalMarket = @(K_to_L,A,L,K_entre,N_entre) K_to_L-(A-K_entre)/max(L-N_entre,1e-4);
% - G/Y is equal to the target, where G = TotalTax-Pensions.
% Note that  Y_corp = F(K_corp,N_corp), where F(.) is Cobb-Douglas
GeneralEqmEqns.govbudget = @(G_to_Y,Pensions,TotalTax,A,Y_entre,K_entre,N_entre,L,alpha) ...
    50*((TotalTax-Pensions)/((A-K_entre)^alpha*(L-N_entre)^(1-alpha)+Y_entre) - G_to_Y);

It turn out that fminsearch was good at setting CapitalMarket close to zero but govbudget was stuck at some positive value. So, instead of changing initial conditions for the GE parameters, we multiplied the equation for govbudget by 50 (!!) and fminsearch converged nicely: both GE conditions got very close to zero.

Is there a way to do this trick with heteroagentoptions.multiGEweights? For example,

heteroagentoptions.multiGEweights=[1,50]

?

2

We tried also fminalgo=5, which does a shooting algorithm.

Problem is that it stopped after a few iterations, leaving the first GE condition equal to 0.09. We would like to force it to be more precise, how can we tighten the converge criteria on the GE equations, if there is any?

Thanks!

3

Yeah.

heteroagentoptions.multiGEweights must be a vector of the same length as the number of general eqm eqns [it is 1-by-length(fieldnames(GeneralEqmEqns)), it is always applied internally and the default is that it is a vector of ones]

The default of toolkit is to minimise the sum-of-squares of the general eqm eqns. The weights are applied after taking squares, before taking sum.

If you want to understand the order of the weights, just look at the GeneralEqmEqns structure that you have created, the weights use the same order [the order is always the order in which you created them]

1 Like
4 
heteroagentoptions.toleranceGEprices=10^(-4); % Accuracy of general eqm prices
heteroagentoptions.toleranceGEcondns=10^(-4); % Accuracy of general eqm eqns

Both the above are their default values.
[should work for any/all fminalgo; tell me if you donā€™t think they are working and I will look]

1 Like
5

Thanks!

The following is unrelated to this topic but refers to the stationary distribution code.

If I set simoptions.verbose = 1, the code returns with this error:

Error using fprintf
Function is not defined for sparse inputs.

Error in StationaryDist_Case1_IterationTan_raw (line 50)
            fprintf('StationaryDist_Case1: after %i iterations the current distance ratio is %8.6f (currdist/tolerance, convergence when reaches 1) \n', counter, currdist/simoptions.tolerance)
            ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
Error in StationaryDist_Case1 (line 261)
        StationaryDistKron=StationaryDist_Case1_IterationTan_raw(StationaryDistKron,PolicyKron,N_d,N_a,N_z,pi_z,simoptions);
                           ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^

Note that I also set simoptions.tolerance = 1e-6.

6

Note sure if they work for fminalgo=5 (shooting)

In HeteroAgentStationaryEqm_Case1 we have

 if isfield(heteroagentoptions,'toleranceGEprices_percent')==0
        heteroagentoptions.toleranceGEprices_percent=10^(-3); % one-tenth of one percent
    end

and then the while loop for the shooting

p_percentchange=Inf;
        while any(p_percentchange>heteroagentoptions.toleranceGEprices_percent) % GeneralEqmConditions>heteroagentoptions.toleranceGEcondns

            p_i=GeneralEqmConditionsFnOpt(p); % using heteroagentoptions.outputGEform=1, so this is a vector (note the transpose)

            GeneralEqmConditionsVec=p_i; % Need later to look at convergence

            % Update prices based on GEstruct following the howtoupdate rules
            p_i=p_i(heteroagentoptions.fminalgo5.permute); % Rearrange GeneralEqmEqns into the order of the relevant prices
            I_makescutoff=(abs(p_i)>heteroagentoptions.updateaccuracycutoff);
            p_i=I_makescutoff.*p_i;

            p_new=(p.*heteroagentoptions.fminalgo5.keepold)+heteroagentoptions.fminalgo5.add.*heteroagentoptions.fminalgo5.factor.*p_i-(1-heteroagentoptions.fminalgo5.add).*heteroagentoptions.fminalgo5.factor.*p_i;

            % Calculate GeneralEqmConditions which measures convergence
            if heteroagentoptions.multiGEcriterion==0
                GeneralEqmConditions=sum(abs(heteroagentoptions.multiGEweights.*GeneralEqmConditionsVec));
            elseif heteroagentoptions.multiGEcriterion==1 %the measure of market clearance is to take the sum of squares of clearance in each market
                GeneralEqmConditions=sqrt(sum(heteroagentoptions.multiGEweights.*(GeneralEqmConditionsVec.^2)));
            end

            % Put new prices into Parameters
            for ii=1:length(GEPriceParamNames)
                Parameters.(GEPriceParamNames{ii})=p_new(ii);
            end

            p_percentchange=max(abs(p_new-p)./abs(p));
            p_percentchange(p==0)=abs(p_new(p==0)); %-p(p==0)); but this is just zero anyway
            % Update p for next iteration
            p=p_new;
        end
        p_eqm_vec=p_new;

From the code above, you can see that the options toleranceGEprices and toleranceGEcondns are never used to control the iterations for the shooting. toleranceGEcondns are never used at all. toleranceGEprices is disregarded because the code sets toleranceGEprices_percent=10^(-3) if the option toleranceGEprices_percent is not present in the option structure.

So as far as I understand, the only stopping criterion for the shooting is when the percentage difference in the GE prices falls below 3 percent. Please correct me if Iā€™m wrong

7 
heteroagentoptions.toleranceGEprices=10^(-4); % Accuracy of general eqm prices
heteroagentoptions.toleranceGEcondns=10^(-4); % Accuracy of general eqm eqns

Should now work for heteroagentoptions.fminalgo=5 (shooting algorithm).

PS. If you have Matlab Optimization Toolbox, heteroagentoptions.fminalgo=8 is lsqnonlin() and it is proper fast while still being very easy to use (no additional setup like for shooting).

1 Like
8

I think the simoptions.verbose=1 error is fixed now (I donā€™t get the error in the code I ran, but I think I understand what it was about and made change to fix that, let me know if it still occurs).

Also, cleaned up the initial guess for stationary dist so it no longer fires up the parallel pool.

1 Like
9

There is a typo that creates an error: heteroagentoptions.heteroagentoptions repeated

Unrecognized field name "heteroagentoptions".

Error in HeteroAgentStationaryEqm_Case1 (line 310)
        while any(p_change>heteroagentoptions.heteroagentoptions.toleranceGEprices) || GeneralEqmConditions>heteroagentoptions.toleranceGEcondns
                           ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
10

Another comment: after fixing the typo, the stopping criteria for the shooting loop are:

while any(p_change>heteroagentoptions.toleranceGEprices) || GeneralEqmConditions>heteroagentoptions.toleranceGEcondns

It would be useful to have an upper bound on the number of iterations as well, otherwise the while loop would go on forever. Note that the code has already the option heteroagentoptions,maxiter, but is never used when fminalgo=5.

11

Fixed typo and added heteroagentoptions.maxiter to fminalo=5 (shooting algo). Thanks.

1 Like
12

We have the following GE condition for capital:

GeneralEqmEqns.CapitalMarket = @(K_to_L,A,L,K_entre,N_entre) K_to_L-(A-K_entre)/(L-N_entre)

The problem is that for some parameter values, N_corp=L-N_entre turns out negative, which does not make sense.

You mentioned to us that it is possible to add a penalty to nudge the routine away from this situation. Could you please show us how you would implement the penalty approach?

Thanks!!

13

Not sure how you would do it with shooting, but for the other fminalgo options you just add in a penalty term

GeneralEqmEqns.CapitalMarket = @(K_to_L,A,L,K_entre,N_entre) K_to_L-(A-K_entre)/(L-N_entre)+(L<N_entre)*10000;

here I just make the penalty constant, +(L<N_entre)*10000, but you could do linear like +(L<N_entre)*(100+10*abs(L-N_entre)) or quadratic +(L<N_entre)*(100+10*abs(L-N_entre)+5*abs(L-N_entre)^2)

PS. You probably need to think if it gets added or subtracted (since goal is to get the general eqm eqn to zero)
PPS. This whole thing is going to get squared, but that is not particularly important.
PPPS. In principle it would be possible to set up the penalty as a whole additional general eqm eqn GeneralEqmEqns.extrapenalty = @(L,N_entre) (L<N_entre)*(100+10*abs(L-N_entre)); but the toolkit may object that there are more general eqm eqns that general eqm params, not sure if it will or not as Iā€™ve never actually tried this.

1 Like
14

Topic: add calibration targets to general equilibrium conditions, and do GE loop and calibration all together

I have to calibrate a parameter of the model so that one model moment matches the corresponding moment in the data. I have two (true) GE conditions (specified earlier in the code, not shown here) and 1 calibration condition. I wrote

Params.data_share_EC = 0.09; %here is the data target
GeneralEqmEqns.target_EC = @(EntreC,Entre,data_share_EC) EntreC/Entre-data_share_EC;
% Specify weights
heteroagentoptions.multiGEweights = [1,1,0.5];

I want to give less weight to the calibration condition. I choose a minimization routine like lsqnonlin which works well in practice.

All is good so far. But what if one of the parameter to calibrate is the persistence or standard deviation of the AR(1) process for the Markov shock z? I donā€™t understand if the toolkit can deal with this case, since z_grid and pi_z are fixed.

1 Like
15

Set up the markov z process using vfoptions.ExogShockFn (and simoptions.ExogShockFn). Idea is that vfoptions.ExogShockFn is a function that takes parameters as inputs, and outputs [z_grid, pi_z]. This is ā€˜parametrizing the exogenous shock processā€™.

Once the shocks are parameterized the rest is standard. Set up the calibration/general eqm that you want, and include the parameter that parametrizes the exogenous shock process as one of the parameters to be determined.

There is also vfoptions.EiidShockFn for i.i.d. e shocks.

Life-Cycle Model A12 (appendix A model 12) in the Intro to Life-Cycle Models shows how to use vfoptions.ExogShockFn (and vfoptions.EiidShockFn).

1 Like
16

Great, thanks! Does ā€œExogShockFnā€ work also for infinite horizon models? I am now calibrating our model with entrepreneursā€¦ (I know infinite horizon does not have iid shocks)

17

Related to my answer above, what I am doing with z is slightly more complicated than the standard setup.

My z variable is a combination of different shocks (labor productivity \varepsilon, entrepreneurial productivity \theta, an iid shock \chi and age).

Here is the code I use to create z and pi_z:

%% Switch to TOOLKIT notation, joint grids to save grid points
% z = (eps,theta,xi,age)
n_z       = [n_eps*n_theta*n_xi+1,1,1,1];
z_grid    = zeros(n_eps*n_theta*n_xi+1,4);
[eps_all,theta_all,xi_all] = ndgrid(eps_grid,theta_grid,xi_grid);
eps_all   = eps_all(:);
theta_all = theta_all(:);
xi_all    = xi_all(:);
% Fill in all rows but the last one
z_grid(1:n_eps*n_theta*n_xi,:) = [eps_all,theta_all,xi_all,age_grid(1)*ones(n_eps*n_theta*n_xi,1)];
% Fill in last row, only for retired (age=R)
z_grid(n_eps*n_theta*n_xi+1,:) = [0,0,0,age_grid(2)]; %when old

% z = (eps,theta,xi,age)
% eps varies first, then theta, then xi
% Note kron in reverse order
p_newborn  = kron(kron(p_newborn_xi,p_unc_theta), p_unc_eps);
pi_young   = kron(kron(pi_xi,pi_theta),  pi_eps);
pi_z       = [(1-p_ret)*pi_young,p_ret*ones(n_eps*n_theta*n_xi,1);
              p_die*p_newborn',  (1-p_die)];

Can I just dump this in vfoptions.ExogShockFn and simoptions.ExogShockFn?

18

Yesterday it didnā€™t, but today ExogShocFn works in infinite horizon :slight_smile:

Yes, just dump all of that into vfoptions.ExogShockFn

1 Like
19

Thanks, very much appreciated!
I was wondering if it is possible to set bounds for the GE variables.

Itā€™s not a priority in any way for me, but it should be possible to do it easily with fminsearch and nonlinear least squares.

fminsearch has the function fminsearchbnd written by John Dā€™Errico which implements Nelder Meade with bounds. I used it extensively in my projects and it works pretty well. Nonlinearlsq allows the user to specify bounds

20

Follow-up on my question above.

I noticed that the toolkit function HeteroAgentStationaryEqm_Case1 now has the option heteroagentoptions.constrainpositive and others to set bounds on parameters during calibration.

My question is actually about the general equilibrium, not calibration/estimation. Suppose I have a model where the GE variables are r and w and I would like to enforce w\geq 0. How do I do that?

The problem is that when I use fminsearch in the toolkit (or similar algorithms chosen with fminalgo), sometimes the code picks a negative w and this generates an error at some point.

Thanks!