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:
- scalar
0or1; - 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.
