Constrain General Eqm Parameters

8

According to my AI, there are a few bugs.

Proofreading notes: parameter-constraint transformations in VFI Toolkit calibration

Focus: parameter transformations used to handle constraints in the calibration routines.

Relevant files:

  • SubCodes/ParameterConstraints/ParameterConstraints_TransformParamsToUnconstrained.m
  • SubCodes/ParameterConstraints/ParameterConstraints_TransformParamsToOriginal.m
  • SubCodes/ParameterConstraints/ParameterConstraints_PType_TransformParamsToUnconstrained.m
  • Estimation/Calibration/CalibrateLifeCycleModel.m
  • Estimation/Calibration/CalibrateLifeCycleModel_PType.m
  • Estimation/Calibration/CalibrateOLGModel_PType.m
  • Estimation/Calibration/CalibrateBIHAModel_PType.m
  • Estimation/Calibration/CalibrateInfHorzAgentModel_PType.m

1. Serious bug: PType constrained parameters are transformed back using the non-PType routine

In PType calibration, the initial vector is transformed using the PType-specific routine:

[calibparamsvec0,caliboptions] = ...
    ParameterConstraints_PType_TransformParamsToUnconstrained( ...
        calibparamsvec0, calibparamsvecindex, CalibParamNames, ...
        nCalibParamsFinder, caliboptions, 1);

This is correct, because PType calibration expands original parameter names into PType-specific blocks.

However, at the end of the PType calibration wrappers, the code transforms back using the non-PType routine:

[calibparamsvec,penalty] = ...
    ParameterConstraints_TransformParamsToOriginal( ...
        calibparamsvec, calibparamsvecindex, CalibParamNames, caliboptions);

This is likely wrong.

The non-PType inverse routine loops over:

for pp = 1:length(CalibParamNames)

But in PType calibration, the relevant number of blocks is:

nCalibParams = size(nCalibParamsFinder,1);

These are not the same object. CalibParamNames contains original parameter names, while nCalibParamsFinder indexes the PType-expanded parameter blocks.

Therefore, if a constrained parameter is PType-specific, some entries may remain in log/logit space instead of being transformed back to the original model scale.

Suggested fix

Add a PType-specific inverse transform:

function [calibparamsvec,penalty] = ...
    ParameterConstraints_PType_TransformParamsToOriginal( ...
        calibparamsvec, calibparamsvecindex, CalibParamNames, ...
        nCalibParamsFinder, caliboptions)

nCalibParams = size(nCalibParamsFinder,1);

penalty = zeros(length(calibparamsvec),1);

for pp = 1:nCalibParams

    idx = calibparamsvecindex(pp)+1:calibparamsvecindex(pp+1);

    if caliboptions.constrainpositive(pp)==1

        temp = calibparamsvec(idx);

        penalty(idx) = abs(temp/50).*(temp < -51);

        calibparamsvec(idx) = exp(temp);

    elseif caliboptions.constrain0to1(pp)==1

        temp = calibparamsvec(idx);

        penalty(idx) = abs(temp/50).*((temp > 51) + (temp < -51));

        calibparamsvec(idx) = 1./(1+exp(-temp));

        calibparamsvec(idx) = calibparamsvec(idx).*(temp > -50);
        calibparamsvec(idx) = calibparamsvec(idx).*(1-(temp > 50)) + (temp > 50);

    end

    if caliboptions.constrainAtoB(pp)==1

        A = caliboptions.constrainAtoBlimits(pp,1);
        B = caliboptions.constrainAtoBlimits(pp,2);

        calibparamsvec(idx) = A + (B-A)*calibparamsvec(idx);

    end

end

if sum(penalty)>0
    penalty = 1/prod(1./penalty(penalty>0));
else
    penalty = 0;
end

end

Then replace, in PType wrappers:

[calibparamsvec,penalty] = ...
    ParameterConstraints_TransformParamsToOriginal( ...
        calibparamsvec, calibparamsvecindex, CalibParamNames, caliboptions);

with:

[calibparamsvec,penalty] = ...
    ParameterConstraints_PType_TransformParamsToOriginal( ...
        calibparamsvec, calibparamsvecindex, CalibParamNames, ...
        nCalibParamsFinder, caliboptions);

Apply this to:

  • CalibrateLifeCycleModel_PType.m
  • CalibrateOLGModel_PType.m
  • CalibrateBIHAModel_PType.m
  • CalibrateInfHorzAgentModel_PType.m

2. Serious bug in CalibrateLifeCycleModel.m: extra usingallstats under fminalgo==8

In the non-lsqnonlin branch, the objective function receives:

..., usingallstats, usinglcp, usingcustomstats, ...

But in the fminalgo==8 branch, it receives:

..., usingallstats, usingallstats, usinglcp, usingcustomstats, ...

There is an extra usingallstats.

Since fminalgo==8 is the default, this is important.

Suggested fix

Change:

..., FnsToEvaluateParamNames, usingallstats, usingallstats, usinglcp, usingcustomstats, targetmomentvec, ...

to:

..., FnsToEvaluateParamNames, usingallstats, usinglcp, usingcustomstats, targetmomentvec, ...

3. Add checks for positive constraints

Current logic:

calibparamsvec(idx) = max(log(calibparamsvec(idx)),-49.99);

This handles p=0, but not p<0. For p<0, log(p) is complex.

Suggested fix

idx = calibparamsvecindex(pp)+1:calibparamsvecindex(pp+1);
p = calibparamsvec(idx);

if any(p < 0)
    error('Initial guess for positive-constrained parameter is negative.');
end

calibparamsvec(idx) = max(log(p), -49.99);

Apply this in:

  • ParameterConstraints_TransformParamsToUnconstrained.m
  • ParameterConstraints_PType_TransformParamsToUnconstrained.m

4. Add checks for 0-to-1 constraints

Current logic:

calibparamsvec(idx) = min(49.99, max(-49.99, ...
    log(calibparamsvec(idx)./(1-calibparamsvec(idx)))));

This handles exactly p=0 and p=1, but not p<0 or p>1.

Suggested fix

idx = calibparamsvecindex(pp)+1:calibparamsvecindex(pp+1);
p = calibparamsvec(idx);

if any(p < 0 | p > 1)
    error('Initial guess for 0-to-1 constrained parameter is outside [0,1].');
end

calibparamsvec(idx) = min(49.99, max(-49.99, log(p./(1-p))));

Apply this in:

  • ParameterConstraints_TransformParamsToUnconstrained.m
  • ParameterConstraints_PType_TransformParamsToUnconstrained.m

5. Add checks for A-to-B constraints

Current logic first maps:

p01 = (p-A)/(B-A);

then applies the 0-to-1 logit transform.

This is valid only if:

  • B > A
  • initial p is in [A,B]

Suggested fix

idx = calibparamsvecindex(pp)+1:calibparamsvecindex(pp+1);

A = caliboptions.constrainAtoBlimits(pp,1);
B = caliboptions.constrainAtoBlimits(pp,2);

if B <= A
    error('constrainAtoB upper bound must be greater than lower bound.');
end

p = calibparamsvec(idx);

if any(p < A | p > B)
    error('Initial guess for A-to-B constrained parameter is outside [A,B].');
end

calibparamsvec(idx) = (p-A)/(B-A);

Apply this in:

  • ParameterConstraints_TransformParamsToUnconstrained.m
  • ParameterConstraints_PType_TransformParamsToUnconstrained.m

6. Penalty is not used in the lsqnonlin residual vector

The inverse transform computes a penalty when transformed parameters move too far beyond the artificial cutoffs.

For scalar objectives, the penalty is used. But for lsqnonlin, the objective returns a vector of residuals with caliboptions.vectoroutput==2. The penalty is not appended to that vector.

Suggested fix

In the vectoroutput==2 branch of each objective function, after constructing Obj, add:

if penalty > 0
    Obj = [Obj; sqrt(penalty)];
end

This is natural because lsqnonlin minimizes the sum of squared residuals.

Apply this to relevant objective functions:

  • CalibrateLifeCycleModel_objectivefn.m
  • CalibrateOLGModel_*_objectivefn.m
  • CalibrateBIHAModel_*_objectivefn.m
  • CalibrateInfHorzAgentModel_objectivefn.m
  • PType objective functions

7. The inverse transformation formula is conceptually fine

The non-PType inverse transform uses:

p = exp(x)

for positive constraints,

p = 1./(1+exp(-x))

for 0-to-1 constraints, and

p = A + (B-A)*p01

for A-to-B constraints.

These formulas are correct.

The main problems are:

  1. PType wrappers call the non-PType inverse transform.
  2. Initial guesses are not checked before applying log or logit.
  3. lsqnonlin does not use the penalty in the residual vector.
  4. CalibrateLifeCycleModel.m has a duplicated usingallstats argument under fminalgo==8.

Summary of recommended changes

Highest priority:

  1. Add ParameterConstraints_PType_TransformParamsToOriginal.m.
  2. Use it in all PType calibration wrappers.
  3. Fix duplicated usingallstats in CalibrateLifeCycleModel.m under fminalgo==8.

Medium priority:

  1. Add explicit initial-guess checks for:
    • positive constraints: parameter must be non-negative
    • 0-to-1 constraints: parameter must be in [0,1]
    • A-to-B constraints: parameter must be in [A,B]
    • A-to-B limits: require B>A

Lower priority:

  1. Add a penalty residual in the vectoroutput==2 / lsqnonlin branch.