Hi Robert,
I recently tried to replicate some papers including calibration. I think you mentioned in one post that we can add calibration targets to GE conditions, but I’m not sure how to do that.
For example, if I want to calibrate beta to target K/H=1.5, what should I put in GE conditions? do I need to treat beta as an additional equilibrium price?
1 Like
In addition to your GE conditions, add one that says
GeneralEqmEqns.KdivHtarget=@(K,H) K/H-1.5;
and then add the parameter being calibrated (beta) to GEPriceParamNames [the list of names of the parameters that are being solved for in the stationary eqm]. Here ‘K’ and ‘H’ will need to have been set up as FnsToEvaluate.
If it is possible to hit all your calibration targets exactly. Done.
If it is not possible to hit all your calibration targets exactly, then it is worth using heteroagentoptions.multiGEweights. By default it is a vector of ones of with length equal to the number of general eqm eqns [heteroagentoptions.multiGEweights=ones(length(fieldnames(GeneralEqmEqns)),1);] and when you solve for general eqm the codes are solving the problem of minimizing the (weighted) sum of squares of the functions in GeneralEqmEqns by choosing the parameters named in GEPriceParamNames.
You can set heteroagentoptions.multiGEweights to any vector of that length. If you look at GeneralEqmEqns you can see the ‘order’ of the names. Set your vector of weights in that same order. You should put larger weights on the actual general equilibrium conditions (as these do have to evaluate to zero) and lower weights on the calibration targets (so that them not evaluating to exactly zero does not ruin your general equilibrium conditions).
PS. If you are calibrating a life-cycle model (no general equilibrium) then there is also the option of using CalibrateLifeCycleModel(). It is not documented anywhere yet, but works in essentially the same way as the GMM estimation commands which are explained in Intro to Life-Cycle Models.
PPS. There are two ways to think about doing calibration (in general eqm models), one is to use joint-optimization (hit calibration targets and general eqm conditions in the same minimization) which is what is happening here. Alternatively, you can do nested-optimization, the inner-optimization would do the general eqm conditions exactly given current parameter, and the outer-optimization would do the calibration targets. Nested will typically be slower, and you get a nice balance between the two by using weights in the joint-optimization approach.
2 Likes
@robertdkirkby If I want to calibrate homeownership rate of retired population, how to modify the FnsToEvaluate?
FnsToEvaluate.oldrate=@() (agej>Jr)*(hprime>0) generates the ratio of house owning of retired population to total population, but I want denominator to be the retired population.
1 Like
Easiest is probably to just first calculate the retired population (as fraction of total population with is mass 1). Because this is exogenous you can just precompute the fraction retired and store it is a parameter, so you have something like
Params.fracretired=0.3;
Then just redefine old rate to be
FnsToEvaluate.oldrate=@() (agej>Jr)*(hprime>0)/fracretired
[following you I omit the input arguments to this]
If you happened to have a model with endogenous retirement you cannot just precompute fracretired as it will (presumably) depend on general eqm and so you would probably want to instead set up both
FnsToEvaluate.oldrate=@() (agej>Jr)*(hprime>0)
and
FnsToEvaluate.fracretired=@() (agej>Jr)
and then in your general eqm condition that you use for calibration you set something like
GeneralEqmEqns.homeownershiprate_retired=@(oldrate,fracretired,target) oldrate/fracretired-target;
[where target is a number that you want to calibrate to, and stored in Params.target]
2 Likes
I have another question related to Chen(2010). The author calibrates beta and theta to match K/Y=1.729, H/Y=1.075. I tried to do it use toolkit and add the following codes:
FnsToEvaluate.ky=@(aprime,hprime,a,h,z,kappaj,r,tau_p,theta,upsilon,phi,alpha,delta_k,delta_o,delta_r,agej,Jr,b,Tr,p) (a-(hprime==0)Chen2010_HousingServicesFn(aprime,hprime,a,h,z,kappaj,r,tau_p,theta,upsilon,phi,alpha,delta_k,delta_o,delta_r,agej,Jr,b,Tr)(1-p))/(((a-(hprime==0)Chen2010_HousingServicesFn(aprime,hprime,a,h,z,kappaj,r,tau_p,theta,upsilon,phi,alpha,delta_k,delta_o,delta_r,agej,Jr,b,Tr)(1-p))^alpha)((kappajz)^(1-alpha)));
FnsToEvaluate.hy=@(aprime,hprime,a,h,z,kappaj,r,tau_p,theta,upsilon,phi,alpha,delta_k,delta_o,delta_r,agej,Jr,b,Tr,p) h/(((a-(hprime==0)Chen2010_HousingServicesFn(aprime,hprime,a,h,z,kappaj,r,tau_p,theta,upsilon,phi,alpha,delta_k,delta_o,delta_r,agej,Jr,b,Tr)(1-p))^alpha)((kappajz)^(1-alpha)));
GEPriceParamNames={‘r’,‘Tr’,‘beta’,‘theta’};
GeneralEqmEqns.capitaloutput=@(ky) ky-1.729;
GeneralEqmEqns.houseoutput=@(hy) hy-1.075;
Params.beta=p_eqm.beta;
Params.theta=p_eqm.theta;
The first two are K/Y = (a-hr*(1-p))/(k^alpha*N^(1-alpha)) and H/Y, and it generates an error:
Error using gpuArray/arrayfun
POWER: needs to return a complex result, but this is not supported for real input X and Y on the GPU. Use POWER(COMPLEX(X),COMPLEX(Y,0)) instead.
Error in EvalFnOnAgentDist_Grid_J (line 96)
Values=arrayfun(FnToEvaluate, PolicyValuesPermute(:,:,:,1),PolicyValuesPermute(:,:,:,2), a_gridvals(:,1),a_gridvals(:,2), z_gridvals_J(1,:,:,1), CellOverAgeOfParamValues{:});
^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
Error in EvalFnOnAgentDist_AggVars_FHorz_Case1 (line 161)
Values=EvalFnOnAgentDist_Grid_J(FnsToEvaluate{ff},ParamCell,PolicyValuesPermute,l_daprime,n_a,n_z,a_gridvals,z_gridvals_J);
^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
Error in chen2010_test (line 246)
AggVars=EvalFnOnAgentDist_AggVars_FHorz_Case1(StationaryDist,Policy, FnsToEvaluate,Params,,n_d,n_a,n_z,N_j,d_grid,a_grid,z_grid,,simoptions);
^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
1 Like
You will want to set up FnsToEvaluate.ky as two separate functions, one for k and one for y. The mean of a ratio is not the ratio of the means. (In the data you have the ratio of the aggregates, but in the model you currently have the mean of the ratio.) FnsToEvaluate commands (AggVars, AllStats, etc.) work by first evaluating that function at every point in the grid, and then calculating statistics (mean, median, variance, Gini, etc.) from them based on the weights in the agent distribution.
That error occurs because there is a negative number being raised to a power, and so the answer would be a complex number [and the GPU does not like this so it throws an error]. It is not immediately obvious to me why you get a negative number, maybe look at what the most recent iteration of the general eqm price parameters is showing and it is perhaps trying a negative value? [I’m guessing the error occurs when you run the general eqm command]
2 Likes
Thanks Robert. So now I use:
FnsToEvaluate.k=@(aprime,hprime,a,h,z,kappaj,r,tau_p,theta,upsilon,phi,alpha,delta_k,delta_o,delta_r,agej,Jr,b,Tr,p) a-(hprime==0)Chen2010_HousingServicesFn(aprime,hprime,a,h,z,kappaj,r,tau_p,theta,upsilon,phi,alpha,delta_k,delta_o,delta_r,agej,Jr,b,Tr)(1-p);
FnsToEvaluate.y=@(aprime,hprime,a,h,z,kappaj,r,tau_p,theta,upsilon,phi,alpha,delta_k,delta_o,delta_r,agej,Jr,b,Tr,p) ((a-(hprime==0)Chen2010_HousingServicesFn(aprime,hprime,a,h,z,kappaj,r,tau_p,theta,upsilon,phi,alpha,delta_k,delta_o,delta_r,agej,Jr,b,Tr)(1-p))^alpha)((kappajz)^(1-alpha));
to represent k and y. Are these correct? The error still exists:
Error using gpuArray/arrayfun
POWER: needs to return a complex result, but this is not supported for real input X and Y on the GPU. Use POWER(COMPLEX(X),COMPLEX(Y,0)) instead.
Error in EvalFnOnAgentDist_Grid_J (line 96)
Values=arrayfun(FnToEvaluate, PolicyValuesPermute(:,:,:,1),PolicyValuesPermute(:,:,:,2), a_gridvals(:,1),a_gridvals(:,2), z_gridvals_J(1,:,:,1), CellOverAgeOfParamValues{:});
^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
Error in EvalFnOnAgentDist_AggVars_FHorz_Case1 (line 161)
Values=EvalFnOnAgentDist_Grid_J(FnsToEvaluate{ff},ParamCell,PolicyValuesPermute,l_daprime,n_a,n_z,a_gridvals,z_gridvals_J);
^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
Error in chen2010_test (line 245)
AggVars=EvalFnOnAgentDist_AggVars_FHorz_Case1(StationaryDist,Policy, FnsToEvaluate,Params,,n_d,n_a,n_z,N_j,d_grid,a_grid,z_grid,,simoptions);
Here is my code. I only changed equilibrium prices and equations.
And
these are variables in workspace when MATLAB stops running. Which one shows the most recent iteration?
1 Like
If you are running HeteroAgentStationaryEqm_Case1_FHorz() with heteroagentoptions.verbose=1 then it prints to the screen (“command window”) the prices, aggvars, and general eqm eqn values at each iteration. You can watch them changing and updating with each iteration as it tries to find the general eqm.
[it won’t be in the workspace as it fails mid-run and doesn’t save the intermediate steps]
1 Like
It is stopping when you call “EvalFnOnAgentDist_AggVars_FHorz_Case1()”. So the issue is one of your FnsToEvaluate is trying to return a complex number (most likely because it does a negative number to a power).
I would comment out your functions to evaluate one-by-one, and run “EvalFnOnAgentDist_AggVars_FHorz_Case1()” each time, that will allow you to figure out which FnsToEvaluate is causing the error. Then look closely at that formula for that one of your FnsToEvaluate and you will probably be able to see what is going wrong.
[E.g., if you normally have
FnsToEvaluate.A=@(aprime,a,z) a;
FnsToEvaluate.Z=@(aprime,a,z,) z;
you could comment one line out, create FnsToEvaluate that just contains one of the two, and see if that works okay in AggVars]
PS. the asset grid in the Chen2010 model can take negative values, so there is a decent chance it is that.
1 Like
You’re right. It turns out the error happens after I add FnsToEvaluate.y, since y=(a-hr(1-p))^alpha*N^(1-alpha). When a is negative or too small, physical capital K becomes negative. Is there any way to bypass this issue?
1 Like
You cannot bypass it in the sense that it is correct, given the FnsToEvaluate.y that you have set.
What you can do is think more carefully about exactly what the variable you want to calculate is, e.g., maybe you should put max(a,0) instead of a. What exactly you want to do about it will depend on exactly what it is you want to calculate.
2 Likes