# Calibrate Bewley models
**URL:** https://discourse.vfitoolkit.com/t/calibrate-bewley-models/667
**Category:** Uncategorized
**Created:** [June 6, 2026, 8:49pm UTC](https://discourse.vfitoolkit.com/t/calibrate-bewley-models/667 "2026-06-06T20:49:13Z")
**Posts on this page:** 16
**Page:** 1
### Author:  [@aledinola](https://discourse.vfitoolkit.com/u/aledinola)
#### Post date: [June 6, 2026, 8:49pm UTC](https://discourse.vfitoolkit.com/t/calibrate-bewley-models/667/1 "2026-06-06T20:49:13Z")
I have some issues with the function to calibrate Bewley models available in the estimation subfolder of the vfi-toolkit.
I was planning to use it to calibrate a model with entrepreneurs I am working on. Is there any example or code that shows how to use it @robertdkirkby? Has it ever been tested?
I will test it on my model but even now my AI found some clear bugs.
# Bugs in `CalibrateBIHAModel.m`
I was looking through `CalibrateBIHAModel.m` and found two issues that look like definite bugs. I have tried to keep this report narrow and only include things that seem unambiguous.
## 1. `calibomitparamsmatrix` assignment fails for non-scalar omitted calibration parameters
In the setup of calibration parameters, the code initializes
```
calibomitparamsmatrix = zeros(1,1);
```
and later, inside the omitted-parameter branch, does
```
tempomitparam = caliboptions.omitcalibparam.(CalibParamNames{pp});
...
calibomitparamsmatrix(1,sum(calibomitparams_counter)) = tempomitparam;
```
The left-hand side is a single scalar location:
```
calibomitparamsmatrix(1,col)
```
but `tempomitparam` can be a vector or array. The surrounding code appears to allow this: omitted entries are identified using `NaN`, and calibrated entries are selected by
```
tempparam = tempparam(isnan(tempomitparam));
```
So this works only when `tempomitparam` is scalar. If `tempomitparam` is non-scalar, MATLAB will error because a vector or array is being assigned into one scalar element.
A minimal reproduction of the issue is:
```
A = zeros(1,1);
x = [1; NaN; 3];
A(1,1) = x;
```
This gives a size-mismatch assignment error.
Later in the function, the code reads the omitted parameter as a column:
```
currparamraw = calibomitparamsmatrix(:,sum(calibomitparams_counter(1:pp)));
currparamraw(isnan(currparamraw)) = calibparamsvec(calibparamsvecindex(pp)+1:calibparamsvecindex(pp+1));
CalibParams.(CalibParamNames{pp}) = currparamraw;
```
So it looks like the intention was to store each omitted-parameter vector as one column of `calibomitparamsmatrix`.
A possible fix would be something like:
```
col = sum(calibomitparams_counter);
calibomitparamsmatrix(1:numel(tempomitparam),col) = tempomitparam(:);
```
Then, when reconstructing the parameter, one may also want to restore the original shape of `tempomitparam`, depending on whether the toolkit wants to preserve row/column orientation.
For example:
```
omitmask = caliboptions.omitcalibparam.(CalibParamNames{pp});
currparamraw = calibomitparamsmatrix(1:numel(omitmask),sum(calibomitparams_counter(1:pp)));
currparamraw = reshape(currparamraw,size(omitmask));
currparamraw(isnan(currparamraw)) = calibparamsvec(calibparamsvecindex(pp)+1:calibparamsvecindex(pp+1));
CalibParams.(CalibParamNames{pp}) = currparamraw;
```
## 2. Non-struct vector `caliboptions.logmoments` branch refers to undefined variables
The function calls
```
[targetmomentvec,usingallstats,usingautocorr,usingcrosssec,usingcustomstats, ...
allstatmomentnames,allstatcummomentsizes,AllStats_whichstats, ...
FnsToEvaluate_AllStats, ...
autocorrmomentnames,autocorrcummomentsizes,AutoCorrStats_whichstats, ...
FnsToEvaluate_AutoCorrStats, ...
crosssecmomentnames,crossseccummomentsizes,CrossSecStats_whichstats, ...
FnsToEvaluate_CrossSecStats, ...
cmsmomentnames,cmscummomentsizes] = SetupTargetMoments_InfHorz(...);
```
So the moment-name and cumulative-size variables available in this function include
```
allstatmomentnames
allstatcummomentsizes
autocorrmomentnames
autocorrcummomentsizes
crosssecmomentnames
crossseccummomentsizes
cmsmomentnames
cmscummomentsizes
```
However, in the non-struct `caliboptions.logmoments` branch, the code refers to variables that do not appear to be defined anywhere in `CalibrateBIHAModel.m`:
```
acsmomentnames
allstatmomentsizes
acsmomentsizes
```
Specifically, the branch contains:
```
if length(caliboptions.logmoments)==(length(acsmomentnames)+length(allstatmomentnames))
temp = caliboptions.logmoments;
caliboptions.logmoments = zeros(length(targetmomentvec),1);
cumsofar = 1;
for mm = 1:length(temp)
if mm <= allstatmomentsizes
caliboptions.logmoments(cumsofar:cumsofar+allstatmomentsizes(mm)) = temp(mm);
cumsofar = cumsofar + allstatmomentsizes(mm);
else
caliboptions.logmoments(cumsofar:cumsofar+acsmomentsizes(mm)) = temp(mm);
cumsofar = cumsofar + acsmomentsizes(mm);
end
end
```
This means that if `caliboptions.logmoments` is non-scalar and has at least one positive entry, MATLAB will try to evaluate
```
length(acsmomentnames)
```
before it ever reaches the later case
```
elseif length(caliboptions.logmoments)==length(targetmomentvec)
```
Since `acsmomentnames` is not defined in the function, this branch errors immediately.
A minimal reproduction of the issue is:
```
caliboptions.logmoments = [1; 0; 0];
if any(caliboptions.logmoments > 0)
if isscalar(caliboptions.logmoments)
% not reached
else
if length(caliboptions.logmoments) == length(acsmomentnames)
% acsmomentnames is undefined
end
end
end
```
A simple conservative fix would be to allow only the two unambiguous non-struct cases:
1. scalar `0` or `1`;
2. a vector already of length `length(targetmomentvec)`.
For example:
```
elseif any(caliboptions.logmoments > 0)
if isscalar(caliboptions.logmoments)
caliboptions.logmoments = ones(length(targetmomentvec),1);
else
if length(caliboptions.logmoments) == length(targetmomentvec)
caliboptions.logmoments = caliboptions.logmoments(:);
else
fprintf('Relevant to following error: length(caliboptions.logmoments)=%i \n', ...
length(caliboptions.logmoments))
fprintf('Relevant to following error: length(targetmomentvec)=%i \n', ...
length(targetmomentvec))
error('You are using caliboptions.logmoments, but it must be scalar or have length equal to the number of target moments.')
end
end
end
```
This avoids the undefined variables and keeps the already-supported vector-of-target-moments case.
---
### Author:  [@aledinola](https://discourse.vfitoolkit.com/u/aledinola)
#### Post date: [June 7, 2026, 11:21am UTC](https://discourse.vfitoolkit.com/t/calibrate-bewley-models/667/2 "2026-06-07T11:21:48Z")
I think there is a more urgent issue: in my Pijoan-Mas test, gridinterplayer now does not work i.e. gives nonsensical results ☹
> **[GitHub - aledinola/PijoanMasTaxes: Aiyagari model with endogenous labor (as in...](https://github.com/aledinola/PijoanMasTaxes)**
>
> Aiyagari model with endogenous labor (as in Pijoan-Mas 2006), with VFI toolkit
---
### Author:  [@aledinola](https://discourse.vfitoolkit.com/u/aledinola)
#### Post date: [June 7, 2026, 11:52am UTC](https://discourse.vfitoolkit.com/t/calibrate-bewley-models/667/3 "2026-06-07T11:52:50Z")
Codex fixed it, sent PR: [Fix grid interpolation policy unpacking with discrete choices by aledinola · Pull Request #104 · vfitoolkit/VFIToolkit-matlab · GitHub](https://github.com/vfitoolkit/VFIToolkit-matlab/pull/104)
@robertdkirkby Let me know if you agree. This change make my Pijoan-Mas replication generate good results.
---
### Author:  [@aledinola](https://discourse.vfitoolkit.com/u/aledinola)
#### Post date: [June 7, 2026, 1:12pm UTC](https://discourse.vfitoolkit.com/t/calibrate-bewley-models/667/4 "2026-06-07T13:12:01Z")
I have a suggestion to improve the outputs of `CalibrateBIHAModel.m`.
Now it displays results on the screen, but for long calibrations (which are the norm in macro!) you would like to save some intermediate results.
Would it be possible to modify `CalibrateBIHAModel.m` so that it saved a log file in txt format? I am using this workflow in my research: every time the objective function improves (i.e. goes down), record an entry in the log file with best parameter values so far and best obj value so far.
**Why this helps**? Because calibrating a model takes usually a lot of time and you rarely let the estimation reoutine complete naturally. As an idea, for a paper of mine where the model takes about 15 seconds to solve, one calibration takes **at least** 2-3 hours to converge. However, a typical minimization routine will improve significantly over maybe the first hour and then it will “stall”, that is, make only tiny improvements. So in practice you want to cut the code after, say, one hour.
Note: I have experienced this with local minimization routines, which are the ones typically used by the toolkit. For simulated annealing is different, since I have seen improving even after several hours.
Glad to share the matlab function I use for this task:
```matlab
function report_estim_progress(obj_smm,Params,calibNames,model_mom,data_mom,moment_contrib,targetNames,write_to_file)
% Print the current estimation state and optionally append the best-so-far one.
%
% INPUTS
% obj_smm : scalar objective value for the current parameter vector
% Params : parameter struct after inserting the current guess
% calibNames : names of calibrated parameters to display
% model_mom : struct of model moments
% data_mom : struct of empirical target moments
% moment_contrib : struct of per-moment objective contributions
% targetNames : names of targeted moments to display
% write_to_file : if true, append this report to estim_best_so_far.txt
%
% OUTPUTS
% none; report is printed to the command window or appended to the txt log
header_params = sprintf('%-20s %12s\n', 'Variable', 'Value');
sep_params = sprintf('%s\n', repmat('-', 1, 35));
header_targets = sprintf('%-20s %-16s %-16s %-16s\n','Moment','Model','Data','Contribution');
sep_targets = sprintf('%s\n', repmat('-', 1, 69));
if ~write_to_file
fprintf('Objective value: %.12f\n\n',obj_smm);
fprintf(header_params);
fprintf(sep_params);
for ii=1:numel(calibNames)
value = Params.(calibNames{ii});
for jj=1:numel(value)
fprintf('%-20s % .16g\n',calibNames{ii},value(jj));
end
end
fprintf('\n');
fprintf(header_targets);
fprintf(sep_targets);
for ii=1:numel(targetNames)
model_vals = model_mom.(targetNames{ii})(:);
data_vals = data_mom.(targetNames{ii})(:);
contrib_vals = moment_contrib.(targetNames{ii})(:);
for jj=1:numel(model_vals)
fprintf('%-20s %-16.8g %-16.8g %-16.8g\n',targetNames{ii},model_vals(jj),data_vals(jj),contrib_vals(jj));
end
end
fprintf('\n');
return
end
res_dir = Params.res_dir;
outfile = fullfile(res_dir,'estim_best_so_far.txt');
fid = fopen(outfile,'a');
if fid==-1
warning('report_estim_progress:OpenFailed', ...
'Could not open best-so-far log file: %s', outfile)
return
end
fprintf(fid,'%s\n',repmat('=',1,68));
fprintf(fid,'Best-so-far estimation result\n');
timestamp_str = char(datetime('now','Format','yyyy-MM-dd HH:mm:ss'));
fprintf(fid,'Timestamp: %s\n',timestamp_str);
fprintf(fid,'Objective value: %.12f\n\n',obj_smm);
fprintf(fid,'%s',header_params);
fprintf(fid,'%s',sep_params);
for ii=1:numel(calibNames)
value = Params.(calibNames{ii});
for jj=1:numel(value)
fprintf(fid,'%-20s % .16g\n',calibNames{ii},value(jj));
end
end
fprintf(fid,'\n');
fprintf(fid,'%s',header_targets);
fprintf(fid,'%s',sep_targets);
for ii=1:numel(targetNames)
model_vals = model_mom.(targetNames{ii})(:);
data_vals = data_mom.(targetNames{ii})(:);
contrib_vals = moment_contrib.(targetNames{ii})(:);
for jj=1:numel(model_vals)
fprintf(fid,'%-20s %-16.8g %-16.8g %-16.8g\n',targetNames{ii},model_vals(jj),data_vals(jj),contrib_vals(jj));
end
end
fprintf(fid,'\n');
fclose(fid);
end %end function
```
---
### Author:  [@aledinola](https://discourse.vfitoolkit.com/u/aledinola)
#### Post date: [June 7, 2026, 1:38pm UTC](https://discourse.vfitoolkit.com/t/calibrate-bewley-models/667/5 "2026-06-07T13:38:34Z")
Another suggestion/request. Now the output `calibsummary` of `CalibrateBIHAModel` is a structure but it contains only one field, `objvalue`. It would be helpful if it contained also the model moments/targets. They can be easily recomputed afterwards but why doing this extra step? `CalibrateBIHAModel` should have them in memory, right?
---
### Author:  [@robertdkirkby](https://discourse.vfitoolkit.com/u/robertdkirkby)
#### Post date: [June 7, 2026, 9:24pm UTC](https://discourse.vfitoolkit.com/t/calibrate-bewley-models/667/6 "2026-06-07T21:24:47Z")
> Would it be possible to modify `CalibrateBIHAModel.m` so that it saved a log file in txt format? I am using this workflow in my research: every time the objective function improves (i.e. goes down), record an entry in the log file with best parameter values so far and best obj value so far.
To be sure I understand, you do not want a list of all attempted parameter vectors with their associated objective fn value. You only want the best one so far?
---
### Author:  [@aledinola](https://discourse.vfitoolkit.com/u/aledinola)
#### Post date: [June 7, 2026, 10:39pm UTC](https://discourse.vfitoolkit.com/t/calibrate-bewley-models/667/7 "2026-06-07T22:39:02Z")
Yes. Whenever the objective function improves, you write current parameter values and objective function value to the txt. That’s the idea. I can share with you the files I am using for a project if you think they might help with the implementation.
This raises the question about where to save the log file. This could be another option hidden in `caliboptions`, maybe?
---
### Author:  [@robertdkirkby](https://discourse.vfitoolkit.com/u/robertdkirkby)
#### Post date: [June 8, 2026, 2:22am UTC](https://discourse.vfitoolkit.com/t/calibrate-bewley-models/667/8 "2026-06-08T02:22:09Z")
> [@aledinola](#):
>
> Codex fixed it, sent PR
The PR was good but turned out the bit above it needs removing too so I have instead just pushed a file that fixes both:
>
and
>
>
> \- RecursiveEqmAggShocks: generalise the Endo/Exo state cascade in
> RecursiveGen…eralEqmWithAggShocks\_InfHorz and MEP\_InfHorz\_Step\_
> MatchExpectations to l\_a in {1..3} and l\_z in {1..5}; thread l\_a/l\_z
> through MatchedExpectationsPath\_InfHorz\_shooting{,\_nod}; gate the
> FnsToEvaluate / AggVarsPath dumps and figure(4)/(5) match heatmaps
> behind recursiveeqmoptions.verbose\>=2; sketch in
> matching\_omitIdiosyncraticExogenousStates.
> \- HeteroAgentStationaryEqm\_Case1\_FHorz\_PType: run PType\_setup\_Parameters
> before ReturnFnParamNamesFn so the per-type Parameters are populated
> when the param-name lookup runs.
> \- jequaloneDist\_PType: treat prod(n\_z)==1 as the no-z case, and switch
> the n\_a x n\_z branch from ndims to all(size(...)==\[n\_a,n\_z\]).
> \- ValueFnIter\_InfHorz\_GridInterpLayer: when N\_d\>0, unkron with
> UnKronPolicyIndexes2\_z so d and aprime are split out, instead of
> UnKronPolicyIndexes1\_z with a \[n\_d,n\_a\] first arg.
>
> Co-Authored-By: Claude Opus 4.7 \
---
### Author:  [@aledinola](https://discourse.vfitoolkit.com/u/aledinola)
#### Post date: [June 8, 2026, 9:19am UTC](https://discourse.vfitoolkit.com/t/calibrate-bewley-models/667/9 "2026-06-08T09:19:54Z")
Great, thanks! I’ve tested on my code and it works well.
---
### Author:  [@aledinola](https://discourse.vfitoolkit.com/u/aledinola)
#### Post date: [June 8, 2026, 12:14pm UTC](https://discourse.vfitoolkit.com/t/calibrate-bewley-models/667/10 "2026-06-08T12:14:47Z")
Another implementation detail: to record only the “best so far” estimates, we need to define a persistent variable. Persistent variables in matlab are more robust than globals (they stay local to the function where they are defined), and indeed they are not discouraged. This is the approach I currently use in my project.
I sent a PR: [Add BIHA calibration progress report by aledinola · Pull Request #105 · vfitoolkit/VFIToolkit-matlab · GitHub](https://github.com/vfitoolkit/VFIToolkit-matlab/pull/105)
Just to explain briefly:
- Added two more fields to `caliboptions`: `progressreport` and `progressreportfilename`; they are set to default values in `CalibrateBIHAModel`.
- In `CalibrateBIHAModel_Joint_objectivefn.m` (line 241) we compute the current scalar objective value, build readable moment names, and call `CalibrateBIHAModel_ProgressReport(...)`. Similar for the nested objective function.
- The helper `CalibrateBIHAModel_ProgressReport` does the job: if `currentObj
### Author:  [@aledinola](https://discourse.vfitoolkit.com/u/aledinola)
#### Post date: [June 8, 2026, 12:59pm UTC](https://discourse.vfitoolkit.com/t/calibrate-bewley-models/667/11 "2026-06-08T12:59:28Z")
I attach the calibration report file that was generated by `CalibrateBIHAModel`.
```auto
CalibrateBIHAModel progress report
Started: 2026-06-08 13:49:12
Only new best objective values are reported.
New best objective
Evaluation: 1
Timestamp: 2026-06-08 13:49:13
Objective: 0.005024649223500906
Parameter values:
beta: 0.945
sigma: 1.458
nu: 2.833
lambda: 0.8560000000000001
theta: 0.64
delta: 0.083
r: 0.037
Model and target moments:
CustomModelStats.corr_h_z: model 0.01217689211305187, target 0.02
CustomModelStats.cv_h: model 0.1998336727608556, target 0.22
CustomModelStats.H: model 0.3564179627812684, target 0.33
CustomModelStats.K_to_Y: model 2.943083637115059, target 3
CustomModelStats.wL_to_Y: model 0.6469328657908955, target 0.64
CustomModelStats.I_to_Y: model 0.2442759418805499, target 0.25
GE price values:
r: 0.037
General equilibrium residuals:
CapitalMarket: -0.002320682790003172
New best objective
Evaluation: 2
Timestamp: 2026-06-08 13:49:15
Objective: 0.0009421654498280088
Parameter values:
beta: 0.9453948262698868
sigma: 1.458
nu: 2.833
lambda: 0.8560000000000001
theta: 0.64
delta: 0.083
r: 0.037
Model and target moments:
CustomModelStats.corr_h_z: model 0.02153280529283889, target 0.02
CustomModelStats.cv_h: model 0.203101726702583, target 0.22
CustomModelStats.H: model 0.3553696228726225, target 0.33
CustomModelStats.K_to_Y: model 3.002929193222047, target 3
CustomModelStats.wL_to_Y: model 0.6396487647203305, target 0.64
CustomModelStats.I_to_Y: model 0.2492431230374299, target 0.25
GE price values:
r: 0.037
General equilibrium residuals:
CapitalMarket: 0.0001170534381692914
New best objective
Evaluation: 17
Timestamp: 2026-06-08 13:49:30
Objective: 7.479035844441271e-06
Parameter values:
beta: 0.9454471654505545
sigma: 1.459316037329529
nu: 2.835059301574856
lambda: 1.053906217134599
theta: 0.6399484477407869
delta: 0.08333337959770608
r: 0.03660534921593452
Model and target moments:
CustomModelStats.corr_h_z: model 0.0181385085798845, target 0.02
CustomModelStats.cv_h: model 0.2195477138132528, target 0.22
CustomModelStats.H: model 0.3316902005670289, target 0.33
CustomModelStats.K_to_Y: model 3.000865229196578, target 3
CustomModelStats.wL_to_Y: model 0.6400800766572171, target 0.64
CustomModelStats.I_to_Y: model 0.2500722412661957, target 0.25
GE price values:
r: 0.03660534921593452
General equilibrium residuals:
CapitalMarket: -4.385112901472005e-05
New best objective
Evaluation: 25
Timestamp: 2026-06-08 13:49:38
Objective: 4.391642810729072e-06
Parameter values:
beta: 0.9455295440806921
sigma: 1.445001310197087
nu: 2.861310434773296
lambda: 1.051660199052065
theta: 0.6400063742581237
delta: 0.08333333338677372
r: 0.03667182219935359
Model and target moments:
CustomModelStats.corr_h_z: model 0.02099065591325617, target 0.02
CustomModelStats.cv_h: model 0.2201251745549516, target 0.22
CustomModelStats.H: model 0.3298618755181313, target 0.33
CustomModelStats.K_to_Y: model 3.001663110261218, target 3
CustomModelStats.wL_to_Y: model 0.6397850577786328, target 0.64
CustomModelStats.I_to_Y: model 0.2501385926821782, target 0.25
GE price values:
r: 0.03667182219935359
General equilibrium residuals:
CapitalMarket: 7.376671332736534e-05
New best objective
Evaluation: 33
Timestamp: 2026-06-08 13:49:46
Objective: 4.021874038689279e-06
Parameter values:
beta: 0.9454937012882687
sigma: 1.448190287031664
nu: 2.861903080378347
lambda: 1.052378374328707
theta: 0.6399969796728858
delta: 0.08333333333338647
r: 0.0366633117278634
Model and target moments:
CustomModelStats.corr_h_z: model 0.01910213166673853, target 0.02
CustomModelStats.cv_h: model 0.2198102954430491, target 0.22
CustomModelStats.H: model 0.3300220327235265, target 0.33
CustomModelStats.K_to_Y: model 2.998375881823546, target 3
CustomModelStats.wL_to_Y: model 0.6402050474599109, target 0.64
CustomModelStats.I_to_Y: model 0.2498646568187881, target 0.25
GE price values:
r: 0.0366633117278634
General equilibrium residuals:
CapitalMarket: -6.936217603098482e-05
New best objective
Evaluation: 41
Timestamp: 2026-06-08 13:49:54
Objective: 3.358294246148939e-07
Parameter values:
beta: 0.9455059712156906
sigma: 1.447473035226541
nu: 2.86059077689366
lambda: 1.052457360977862
theta: 0.6400000717291408
delta: 0.08333295852492338
r: 0.03666753091632199
Model and target moments:
CustomModelStats.corr_h_z: model 0.02000671659048116, target 0.02
CustomModelStats.cv_h: model 0.2200055248474263, target 0.22
CustomModelStats.H: model 0.3299763599838959, target 0.33
CustomModelStats.K_to_Y: model 3.000530858044369, target 3
CustomModelStats.wL_to_Y: model 0.639934837687285, target 0.64
CustomModelStats.I_to_Y: model 0.2500431135461642, target 0.25
GE price values:
r: 0.03666753091632199
General equilibrium residuals:
CapitalMarket: 2.174391169820766e-05
New best objective
Evaluation: 49
Timestamp: 2026-06-08 13:50:02
Objective: 5.590842512119035e-08
Parameter values:
beta: 0.9455030919387456
sigma: 1.447402839215763
nu: 2.859677233102068
lambda: 1.052655449924112
theta: 0.6399998124757905
delta: 0.0833333329442448
r: 0.03666647808097931
Model and target moments:
CustomModelStats.corr_h_z: model 0.01999221555249907, target 0.02
CustomModelStats.cv_h: model 0.2200223687946467, target 0.22
CustomModelStats.H: model 0.3299941533407194, target 0.33
CustomModelStats.K_to_Y: model 2.999784537947013, target 3
CustomModelStats.wL_to_Y: model 0.6400264238667074, target 0.64
CustomModelStats.I_to_Y: model 0.249982043661736, target 0.25
GE price values:
r: 0.03666647808097931
General equilibrium residuals:
CapitalMarket: -8.870588484699571e-06
```
---
### Author:  [@aledinola](https://discourse.vfitoolkit.com/u/aledinola)
#### Post date: [June 12, 2026, 4:02pm UTC](https://discourse.vfitoolkit.com/t/calibrate-bewley-models/667/12 "2026-06-12T16:02:33Z")
Hi @robertdkirkby, what do you think about the log file that records estimation/calibration progress? I find it very useful for long calibration runs.
---
### Author:  [@aledinola](https://discourse.vfitoolkit.com/u/aledinola)
#### Post date: [June 21, 2026, 8:00pm UTC](https://discourse.vfitoolkit.com/t/calibrate-bewley-models/667/13 "2026-06-21T20:00:00Z")
Another question regarding calibration.
What happens if any of the model moments in the calibration targets is nonfinite (inf or nan)? If this happens, the risk is that the inf or nan value propagates and the scalar objective value becomes inf or nan.
Looking at the function that computes the objective (see [VFIToolkit-matlab/Estimation/ObjectiveFn/CalibrateBIHAModel\_Joint\_objectivefn.m at master · vfitoolkit/VFIToolkit-matlab · GitHub](https://github.com/vfitoolkit/VFIToolkit-matlab/blob/master/Estimation/ObjectiveFn/CalibrateBIHAModel_Joint_objectivefn.m)),
It seems that when the distance is the sum of squares, the nan is omitted:
```auto
Obj1=sum(caliboptions.weights.* ...
(currentmomentvec(actualtarget)-targetmomentvec(actualtarget)).^2,'omitnan');
```
since the `sum` function is called with the option `omitnan`. This is an acceptable choice but is not the only one: if one set of parameters delivers a moment that is “strange”, maybe the code should assign a large penalty, instead of ignoring the strange moment. In this way the code is pushed away from the “bad” set of parameters.
However, `Inf` is not omitted by ‘omitnan’. If the residual is `Inf` and the weight is positive, the objective becomes `Inf`. This might crash the code: it depends on how the matlab solver behaves with a `Inf` value.
My suggestion: assign a large but finite penalty over allowing `NaN` or `Inf` to reach the optimizer.
---
### Author:  [@robertdkirkby](https://discourse.vfitoolkit.com/u/robertdkirkby)
#### Post date: [June 27, 2026, 3:14pm UTC](https://discourse.vfitoolkit.com/t/calibrate-bewley-models/667/14 "2026-06-27T15:14:42Z")
> [@aledinola](#):
>
> What happens if any of the model moments in the calibration targets is nonfinite (inf or nan)? If this happens, the risk is that the inf or nan value propagates and the scalar objective value becomes inf or nan.
‘NaN’ is used by VFI Toolkit as a mask for calibration/estimation targets. So say for example you have a model with 5 periods, the first three are working age and the last two are retirement. You want to estimate the model to target mean hours worked during working ages. You do this by creating
`TargetMoments.AgeConditionalStats.HoursWorked.Mean=[2,3,3,nan, nan];`
(and create `FnsToEvaluate.HoursWorked`)
The toolkit understands that it should calculate the age-conditional mean of hours worked, and choose parameters to get this close to the target value. The two nan are interpreted as saying that these two ‘ages’ should be ignored during the calibration.
This behaviour of nan as a mask is of course what happens if you choose to put nan into the TargetMoments. It is not what happens if the model returns nan for the value of those moments (given current parameters).
I don’t know exactly how the different optimizers handle `Inf`, but I have not had an issue with this in the past.
Your suggestion over assigning a large but finite penalty over allowing NaN or Inf is reasonable. Before overwriting Inf I would want to see an example where passing Inf causes problems (I would want to see that the optimizers don’t handle Inf already). In the case of NaN it might make more sense to simply error in the case that the model moment is NaN, as this is likely related to parameter values and simply adding parameter constraints would be a more sensible approach that obscuring it from the user.
PS. I am still planning to do the recording progress of calibration/estimation, just didn’t have time in the day or two before I left so will get done when I get back early July.
---
### Author:  [@aledinola](https://discourse.vfitoolkit.com/u/aledinola)
#### Post date: [June 27, 2026, 6:05pm UTC](https://discourse.vfitoolkit.com/t/calibrate-bewley-models/667/15 "2026-06-27T18:05:05Z")
Thanks for your answer.
In a model with entrepreneurs, I have encountered several cases where some simulated moments become nonfinite. The issue is that, for some parameter values, the model can converge to a solution with no entrepreneurs, meaning that all agents choose to be workers. In that case, moments with the mass of entrepreneurs in the denominator can become `Inf`. For example, in MATLAB, `34.0/0` returns `Inf`.
There is a related issue when both the numerator and denominator of a moment are zero. For example, if there are no entrepreneurs, the numerator of some entrepreneur-specific moments may also be zero. In that case, MATLAB returns `NaN`, since `0/0` is `NaN`.
These nonfinite moments can then propagate into the calibration objective, turning the objective value into `Inf` or `NaN`. If the objective is `Inf`, this is not necessarily a major problem: it effectively tells the optimizer to move away from that region of the parameter space. However, if the objective is `NaN`, many routines cannot rank or compare the objective value, and the calibration may fail or stop.
If I understand correctly, you are saying that the toolkit should not hide this problem, because encountering such cases should encourage the user to impose tighter calibration bounds, @robertdkirkby. I agree with that principle. However, in practice it can be difficult to know the relevant boundaries of the parameter space in advance. In my case, a solution with zero entrepreneurs is not mathematically invalid; it is simply a very poor fit for calibration purposes.
For this reason, my experience is that it is useful to prevent the routine from failing completely when a nonfinite moment is produced. Ideally, the code would record the parameter vector and return a large penalty value, rather than allowing a `NaN` to stop the calibration. This is especially important if there is no log file recording the progress made up to that point (so thanks again for considering this!)
**Update**
I looked into the documentation of fminsearch (default fminalgo in the toolkit) and found that with NaN, `fminsearch` may continue, but the result is **not reliable**. MathWorks explicitly warns that solvers expect real objective values, and `NaN` or complex values can cause “unexpected results.”
Inf is bad but can still be ordered relative to finite values. NaN is worse because comparisons involving NaN are false, so it can interfere with the Nelder-Mead logic inside fminsearch.
---
### Author:  [@aledinola](https://discourse.vfitoolkit.com/u/aledinola)
#### Post date: [June 27, 2026, 6:07pm UTC](https://discourse.vfitoolkit.com/t/calibrate-bewley-models/667/16 "2026-06-27T18:07:16Z")
> [@robertdkirkby](#):
>
> PS. I am still planning to do the recording progress of calibration/estimation, just didn’t have time in the day or two before I left so will get done when I get back early July.
I actually implemented this in a branch of the vfi toolkit repo. The branch is called ‘ale’. You can have a look at my implementation if it might help. There is a non-trivial issue with persistent variables. I think I handled it correctly but be aware 🙂