Calibration and Estimation with global optimizers

1

Sometimes (more than sometimes), the estimation of a structural model is challenging, especially if there are many parameters to be calibrated and/or the model has important non-linearities. The objective function can be non-differentiable which makes optimizers like fmincon not well behaved.

The toolkit default uses fminsearch, which is a derivative-free (but still local) method based on Nelder-Meade and is a good compromise between speed and robustness. For difficult estimations, however, one needs a global optimizer. The toolkit has cma-es (covariance matrix adaptation). Matlab provides several built-in global optimizers such as:

and many others. In my esperience I have had good results with cma-es and simulated annealing as well, but I’d like to know what other users think about this.

(I am not suggesting to add any of the matlab routines from the global optim toolbox to the toolkit as of now, just brainstorm).

Of course, the main ingredient of structural estimation is… patience :slight_smile:

1 Like
2

One of my favorite TEDx talks came from a woodsman in North Carolina who said

We exploit the weakness of wood when we work it. We exploit the strength of wood when we use it.

He demonstrates with his axe just how easily wood can be worked into whatever shape (to make a chair) and then how strong the chair can be when loaded. I think that a good toolkit, and a good tool user, knows how to align the tools with the grain of the problem so that the problem yields easily, and that the final product is robust.

When I see gradient descent solve a problem in seconds that takes fminsearch hours and hours to solve, I know that patience is not the best tool…

1 Like
3

I find CMA-ES is what I turn to when things get hard, as while slow it tends to be realiable when all the other algorithms fail.

But it is not the first thing I turn to, I typically start with lsqnonlin() [fminalgo=8 in VFI Toolkit] as when it works it is fast (it depends on derivatives). Sometimes this works, especially if the initial guess is decent, but sometimes it fails and then I turn to CMA-ES.

My other main trick in practice is to start with a smaller sub-model. E.g., I might turn off the shocks and calibrate that for half the parameters using half the targets, then turn the shocks back on and use the submodel parameters as my initial guess. One of the nice things about VFI Toolkit is this kind of thing is genuinely easy, you just set n_z=1 and done, similarly I might exogenise labor by setting n_d=1. Obviously the submodels run way faster, but because you might only calibrate half the parameters and targets, the optimization problem is also much easier.

I have often thought about surrogate optimization as a way to calibrate/estimate larger models. I’m yet to really see a paper that shows that it is or isn’t a good idea. I’d be very interested to see something systematic about using it for heterogeneous agent models.

Slightly off-topic, I was once watching a talk by Michael Keane in which he mentioned his 2004 paper with Susumu Imai, Intertemporal Labor Supply and Human Capital Accumulation. I’ve met Susumu and he tells me the paper took a decade to write and publish. Keane in this talk was saying that Susumu set the maximum likelihood estimation of the life-cycle model off to run, then sat down and spent the month reading Buddenbrooks by Thomas Mann while it converged. [Personally I don’t have the patience for working with models that take this long, but I do like the story.]

PS. The other alternative is of course Bayesian estimation. This side steps the optimization problem entirely, instead you just have to simulate MCMC for long enough. This is not something the toolkit will do for you though (if anyone out there wants to be in charge of Bayesian estimation in VFI Toolkit, please get in touch with me, I’d love to see it but I don’t want to have to learn all the knowledge required to do Bayesian estimation, I have plenty to do with implementing toolkit features already).

2 Likes
4

Yes, but it depends on the properties of the function you have to minimize, no? Gradient descent is great but the function must be differentiable (at least in the basic version of gradient descent I know). Typically, the objective function in a GMM or SMM is not smooth.

1 Like
5

Half guessing, but I don’t think the problem in calibration/estimation of these models is so much smoothness as it is local vs global optima, my impression is many models have large numbers of local optima you can get stuck in.

2 Likes
6

Came across this article by Chen, Didisheim & Scheidegger (2026) who do surrogate estimation of option-pricing models (in Finance), constructing the surrogate using Deep Neural Nets, and argue that it works well including for GMM estimation. [Although like all these Swiss papers they have insane compute facilities, like 10-100x bigger than what everyone else has.]

I looked at Matlab surrogateopt() and it says the surrogate is constructed by first evaluating original model at lots of points (all surrogate approaches share this first step) and then “Create a surrogate model of the objective function by interpolating a radial basis function through all of the random trial points”. Not sure how well this works vs DNNs. Clearly using DNNs for surrogates would require a lot of coding (although Matlab does have good DNN functions so they are fairly easy to build).

2 Likes
7

Hi Robert, do you have an example of the estimation of the parameters of a structural model with the toolkit, preferably from a published paper?

8

The Calibrate OLG command is used in Kirkby, Ngoc & Yao (2026).

The Calibrate Life Cycle model (with permanent types) command is used in example of Huggett, Ventura & Yaron (2006).

The GMM Estimation of Life-Cycle model command is used in replication of Gourinchas & Parker (2002).

You can also find a bunch of examples showing different features of the GMM Estimation of Life-Cycle model commands, as the Life-Cycle Models 46-50 in the Intro to Life-Cycle Models.

Note that other than the calculating of confidence intervals, standard errors, and other measures of parameter uncertainty (and hence needing the covariance matrix of data moments as an input), the GMM estimation and calibration commands are very similar.

There are not yet any examples of the calibration commands for infinite horizon models, hopefully they will be available next month. But the way they function is similar.

3 Likes
9

The estimation of the parameters and the general equilibrium are done in the same command? Is it similar to a nested fixed point?

I had a look at the replication of CastanedaDiazGimenezRiosRull2003.m and I was not able to figure out.

10

Short Answer:
BIHA and OLG calibration commands are (by default) based on joint-optimization of the general eqm and calibration targets.

Long Answer:
Let \theta be a vector of parameters. We can split \theta=[\theta_1, \theta_2] into \theta_1 which is a vector parameters to be determined in general eqm (GeneralEqmParamNames in most codes gives their names) and \theta_2 which is a vector of parameters to be calibrated (CalibParamNames is in most codes gives their names).

We will first describe solving stationary general eqm and calibration separately, and then how to do them together which is what we are doing when using the OLG and BIHA calibration commands.

Let \mathcal{M} be our model (OLG, BIHA, etc.).

First, how to solve stationary general eqm. We want the model to solve n general equilibrium equations, GE_1(\mathcal{M}(\theta)), GE_2(\mathcal{M}(\theta)),…,GE_n(\mathcal{M}(\theta)) [GeneralEqmEqns in most codes, they will evaluate to zero in stationary general eqm]. So solving stationary general eqm (using HeteroAgentStationaryEqm commands) involves solving,
\begin{equation} \min_{\theta_1 \in \Theta_1} \quad \sum_{g=1}^n w_g GE_g(\mathcal{M}(\theta_1))^2 \end{equation}
so minimizing the sum-of-squares of the general eqm eqns. By default the weights, w_i are all set to one.
[By default VFI Toolkit solves stationary general eqm by minimizing the sum-of-squares of the GeneralEqmEqns, you can instead set heteroagentoptions.multiGEcriterion=0 (default is =1) to get sum-of-absolute-values if you want. I don’t know any situation you would actually want to do this, but the option is there.]
[If you want to change the weights, you can set heteroagentoptions.multiGEweights to be a vector of length n, by default it is internally set to heteroagentoptions.multiGEweights=ones(1,length(fieldnames(GeneralEqmEqns)))]

Now we turn to calibration. At first, pretend that we are just calibrating a model without having to worry about stationary general eqm, e.g., calibrating a life-cycle model.

Calibrating the model involves some calibration targets, e.g., the mean of assets, and the codes will aim to get model moments/statistics to hit these targets. Say our targets are A_1,A_2,…,A_m. By the way we name TargetStats in the code the toolkit understands how to compute these from the model as functions a_1(\mathcal{M}(\theta)), a_2(\mathcal{M}(\theta)),…, a_m(\mathcal{M}(\theta)). VFI Toolkit performs the calibration by solving,
\begin{equation} \min_{\theta_2 \in \Theta_2} \quad \sum_{c=1}^m w_c (A_c-a_c(\mathcal{M}(\theta_2)))^2 \end{equation}
You can control the weights w_c using caliboptions.weights. You can change from default of sum-of-squares (caliboptions.metric=‘sum_squared’) to instead using log-ratio (=‘sum_logratiosquared’). You can use caliboptions.logmoments to target the log of moments of the model (you can name the TargetStats, and only log some but not others).

So far so good. Now how to do both calibration and stationary general equilibrium at the same time. There are two conceptual approaches.

The first is ‘joint-optimization’, which is simply to solve,
\begin{equation} \min_{[\theta_1,\theta_2] \in \Theta} \quad w_{rel} \sum_{g=1}^n w_g GE_g(\mathcal{M}(\theta))^2 + \sum_{c=1}^m w_c (A_c- a_c(\mathcal{M}(\theta)))^2 \end{equation}
so we are just doing the two problems jointly. Notice there is w_{rel} equals 10 by default and you can set using caliboptions.relativeGEweight.

The other is ‘nested-optimization’, in this we solve,
\begin{equation} \min_{\theta_2 \in \Theta_2} \quad \sum_{c=1}^m w_c (A_c- a_c(\mathcal{M}(\theta_2)))^2 + F_{GE}(\theta_2) \end{equation}
subject to
\begin{equation} F_{GE}(\theta_2) = \min_{\theta_1 \in \Theta_1} \quad \sum_{g=1}^n w_g GE_g(\mathcal{M}([\theta_1, \theta_2]))^2 \end{equation}
notice that this means one ‘outer’ optimization (for calibration) and one ‘inner’ optimization (for general eqm). So every time we try one a value of the parameter vector \theta_2 for calibration we have to solve the entire stationary general eqm problem.

By default VFI Toolkit does joint-optimization (caliboptions.jointoptimization=1), but you can set caliboptions.jointoptimization=0 to get nested-optimization if you want.

In practice, there is no real reason to use nested-optimization, because it is painfully slow. So much so that while the joint-optimization codes have been tested on numerous examples I haven’t really bothered to actually test the joint-optimization codes, just wrote them [I don’t have the spare compute to be bothered waiting for the nested-optimization codes to run].

After you calibrate a model with stationary general eqm (like OLG or BIHA) you should run stationary general eqm command once before using the model for anything. The calibration command (because it is joint-optimization) will get the general eqm eqns close-ish to zero, but balance this against the calibration targets. But your model requires the general eqm eqns to hold, not be balanced off against things. So you need to run the stationary general eqm command to get the general eqm to hold (this will come at some minor cost to your calibration targets, but typically tiny enough to just ignore). Note that the use of the weights (both in heteroagentoptions, caliboptions, and the relativeGEweight) can be important for getting the joint-optimization to choose parameters that make the general eqm eqns near zero (so running the stationary eqm command once more after calibration doesn’t much ruin your calibration targets); putting enough weight on general eqm eqns that optimizer makes them near-zero, but not enough that it ignores calibration targets.

PPS. In the FHorz calibration commands you can target anything in AllStats or AgeConditionalStats. In the InfHorz calibration commands you can target anything in AllStats, AutoCorrTransProbs, or CrossSectionCovarCorr.

PPPS. You can target some very complex model statistics using CustomModelStats (or equally use them in general eqm eqns). Just requires you to hand-write the code that calculates the complicated model statistics.

PPPPS. Obviously the model solution will (normally) depend on all the parameters in the model. I above omit from \mathcal{M}(\cdot) all the parameters that are not changing as part of solving the calibration/equilibrium. I suspect this is obvious, but feel I should mention it anyway.

1 Like
11

The CastanedaDiazGimenezRiosRull2003.m lacks a functional subroutine for calibration.

Undefined function 'CalibrateFromModelTargetsFn' for input arguments of
type 'function_handle'.

Error in CastanedaDiazGimenezRiosRull2003 (line 526)
[Params1,fval,counteval,exitflag]=CalibrateFromModelTargetsFn(Params,
ParamNamesToEstimate, EstimationTargets, ModelTargetsFn, estimationoptions,
GEPriceParamNames, GeneralEqmTargetNames);
^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^

PS: The 135 line of the EvalFnOnAgentDist_AutoCorrTransProbs_InfHorz did not change the name of the routine: PolicyInd2Val_Case1.

Undefined function 'PolicyInd2Val_Case1' for input arguments of
type 'gpuArray'.

Error in EvalFnOnAgentDist_AutoCorrTransProbs_InfHorz (line 135)
PolicyValues=PolicyInd2Val_Case1(Policy,n_d,n_a,n_z,d_grid,a_grid,simoptions);
^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
Error in CastanedaDiazGimenezRiosRull2003 (line 314)
CorrTransProbs=EvalFnOnAgentDist_AutoCorrTransProbs_InfHorz(StationaryDist,Policy,FnsToEvaluate3,Params,[],n_d,n_a,n_z,d_grid,a_grid,z_grid,pi_z,simoptions);
^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
2 Likes
12

Fixed. Thanks, turns out there were a handful left, searched the repo and cleaned them up.

1 Like
13

When I was doing the CDGRR2003 replication I intended to replicate the calibration. But it turn out the weights and objective function used have been lost to the sands of time. Hence I just abandoned the attempt and the code there is just what was done at the time.

That said, nowadays there is the Calibrate BIHA model command that can calibrate this model. I might go rewrite it just so people can see how you would go about calibrating the model.

PS. Just realised that CDGRR2003 is what @jake88 looked at, so no surprises he couldn’t see how it worked as it doesn’t work :rofl: oops. Sorry Jake, I didn’t put two and two together when you said that is the one you looked at.

14

Just pushed an update of the Castaneda, Diaz-Gimenez & Rios-Rull (2003) replication.

Moved it over to using vfoptions.ExogShockFn (that it wasn’t already shows how old it was, and how much better and easier the toolkit is now).

Also made small improvements to the codes that simulate the two painful moments.

Cleaned up the bottom to use CalibrateBIHAModel command for the calibration. Likely works, I only let it run a few iterations and it was going seemingly okay. Doing it properly probably means improving the weights it uses during calibration. But the thing runs now.

Is actually a nice illustration of how easy the toolkit makes model calibration. It’s calibrating 25 parameters for 28 targets, and yet the whole setup is pretty simple (albeit with a heavy dependence on CustomModelStats).

[99% of the runtime from each calibration iteration is in calculating the two awkward statistics. Solving the model and all other 26 statistics together make up the remaining 1%. So it is kind of silly. But this is just a reflection of the nature of two of the statistics.]

2 Likes
15

I ran the script and received the following error messages.

Error using ValueFnIter_InfHorz (line 211)
Problem with pi_z in ValueFnIter_Case1: min(min(pi_z))<0 \n

Error in CalibrateBIHAModel_Joint_objectivefn (line 49)
    [V, Policy]=ValueFnIter_InfHorz(n_d,n_a,n_z,d_grid, a_grid, z_gridvals, pi_z, ReturnFn, Parameters, DiscountFactorParamNames,[], vfoptions);
                ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
Error in CalibrateBIHAModel>@(calibparamsvec)CalibrateBIHAModel_Joint_objectivefn(calibparamsvec,CalibParamNames,n_d,n_a,n_z,d_grid,a_grid,z_gridvals,pi_z,ReturnFn,Parameters,DiscountFactorParamNames,GEPriceParamNames,ParametrizeParamsFn,FnsToEvaluate,GeneralEqmEqnsCell,GEeqnNames,GeneralEqmEqnParamNames,usingallstats,usingautocorr,usingcrosssec,usingcustomstats,targetmomentvec,allstatmomentnames,autocorrmomentnames,crosssecmomentnames,cmsmomentnames,allstatcummomentsizes,autocorrcummomentsizes,crossseccummomentsizes,cmscummomentsizes,AllStats_whichstats,AutoCorrStats_whichstats,CrossSecStats_whichstats,FnsToEvaluate_AllStats,FnsToEvaluate_AutoCorrStats,FnsToEvaluate_CrossSecStats,nCalibParams,calibparamsvecindex,calibomitparams_counter,calibomitparamsmatrix,caliboptions,heteroagentoptions,vfoptions,simoptions) (line 361)
        CalibrationObjectiveFn=@(calibparamsvec) CalibrateBIHAModel_Joint_objectivefn(calibparamsvec, CalibParamNames,n_d,n_a,n_z,d_grid, a_grid, z_gridvals, pi_z, ReturnFn, Parameters, DiscountFactorParamNames, GEPriceParamNames, ParametrizeParamsFn, FnsToEvaluate, GeneralEqmEqnsCell, GEeqnNames, GeneralEqmEqnParamNames, usingallstats,usingautocorr,usingcrosssec,usingcustomstats, targetmomentvec, allstatmomentnames,autocorrmomentnames,crosssecmomentnames,cmsmomentnames, allstatcummomentsizes,autocorrcummomentsizes,crossseccummomentsizes,cmscummomentsizes, AllStats_whichstats,AutoCorrStats_whichstats,CrossSecStats_whichstats, FnsToEvaluate_AllStats, FnsToEvaluate_AutoCorrStats, FnsToEvaluate_CrossSecStats, nCalibParams, calibparamsvecindex, calibomitparams_counter, calibomitparamsmatrix, caliboptions, heteroagentoptions, vfoptions,simoptions);
                                                 ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
Error in snls (line 340)
            newfvec = feval(funfcn{3},xcurr,varargin{:});
                      ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
Error in lsqncommon (line 178)
            snls(funfcn,xC,lb,ub,flags.verbosity,options,defaultopt,F,Jac,caller, ...
            ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
Error in lsqnonlin (line 264)
    lsqncommon(funfcn,xCurrent,lb,ub,options,defaultopt,caller,...
    ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
Error in CalibrateBIHAModel (line 421)
    [calibparamsvec,calibobjvalue]=lsqnonlin(CalibrationObjectiveFn,calibparamsvec0,[],[],[],[],[],[],[],minoptions);
                                   ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
Error in CastanedaDiazGimenezRiosRull2003 (line 577)
[CalibParams,calibsummary]=CalibrateBIHAModel(CalibParamNames,TargetMoments,n_d,n_a,n_z,d_grid, a_grid, z_grid, pi_z, ReturnFn, Params, DiscountFactorParamNames, [], GEPriceParamNames, FnsToEvaluate, GeneralEqmEqns, heteroagentoptions, caliboptions, vfoptions,simoptions);
                           ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
1 Like