OccBin Toolbox  V2.0

The model and its constraint

The model is a standard RBC model with a capital adjustment cost (PSI) and a Lagrange multiplier lambdak on the investment floor. In the reference regime (constraint not binding), the multiplier is zero: lambdak = 0. In the alternative regime (constraint binding), investment is pegged to its lower bound: i = log(PHII · ISS).

The occasionally binding constraint is I_t ≥ φ · I_ss, or in log-deviation form: ĩ_t ≥ log(φ). The constraint is violated — and regime switches — when the linear solution produces investment below this floor.

Variables in both .mod files:

The two .mod files

The reference file dynrbc.mod sets lambdak = 0 (equation 4 below). The alternative file dynrbcirr_i.mod is an exact replica except that equation 4 is replaced by i = log(PHII · ziss), enforcing the investment floor. All other equations and the parameter list are identical.

// dynrbc.mod — reference regime (constraint not binding)
var a, c, i, k, kprev, lambdak;
varexo erra;
parameters ALPHA, DELTAK, BETA, GAMMAC, RHOA, PHII, PHIK, PSI, PSINEG, ISS, KSS;

model;
// 1. Euler equation (capital)
-exp(c)^(-GAMMAC)*(1+2*PSI*(exp(k)/exp(k(-1))-1)/exp(k(-1)))
+ BETA*exp(c(1))^(-GAMMAC)*((1-DELTAK)
  -2*PSI*(exp(k(1))/exp(k)-1)*(-exp(k(1))/exp(k)^2)
  +ALPHA*exp(a(1))*exp(k)^(ALPHA-1))
  = -lambdak + BETA*(1-DELTAK+PHIK)*lambdak(1);

// 2. Resource constraint
exp(c)+exp(k)-(1-DELTAK)*exp(k(-1))+PSI*(exp(k)/exp(k(-1))-1)^2
  = exp(a)*exp(k(-1))^ALPHA;

// 3. Investment definition
exp(i) = exp(k)-(1-DELTAK)*exp(k(-1));

// 4. Reference regime: multiplier is zero
lambdak = 0;

// 5. TFP process
a = RHOA*a(-1) + erra;

kprev = k(-1);
end;

// dynrbcirr_i.mod — alternative regime (investment at floor)
// Identical to dynrbc.mod except equation 4:

// 4. Alternative regime: investment is pinned to its lower bound
i = log(PHII*ziss);

Parameter values

Parameter values are set in paramfile_dynrbc.m, loaded by the Dynare steady-state file at parse time.

BETA   = 0.96;    % discount factor
ALPHA  = 0.33;    % capital share
DELTAK = 0.10;    % depreciation rate
GAMMAC = 2;       % inverse elasticity of substitution
RHOA   = 0.9;     % TFP persistence
PHII   = 0.975;   % investment floor as fraction of I_ss
PSI    = 0;       % capital adjustment cost (zero here)

The toolkit computes steady-state capital and investment from these primitives: K_ss = ((1/β − 1 + δ) / α)^1/(α−1),  I_ss = δ K_ss.

Setting up and calling the solver

The script runsim_irrcap.m specifies the shock sequence and calls solve_one_constraint. Here, a positive 2% TFP shock hits at period 10 and a symmetric negative shock at period 30, over a 100-period horizon.

clear
global M_ oo_

modnam     = 'dynrbc';        % reference regime .mod file
modnamstar = 'dynrbcirr_i';   % alternative regime .mod file

% Constraint string (in deviation from steady state)
constraint       = 'i<log(PHII)';   % switch to alt. regime when i hits floor
constraint_relax = 'lambdak<0';     % return to ref. regime when multiplier goes negative

irfshock = char('erra');       % TFP innovation

shockssequence = zeros(60,1);
shockssequence(10) =  0.02;      % positive shock at period 10
shockssequence(30) = -0.02;      % negative shock at period 30

nperiods = 100;
maxiter  = 10;

[zdatalinear, zdatapiecewise, zdatass, oobase_, Mbase_] = ...
  solve_one_constraint(modnam, modnamstar, ...
    constraint, constraint_relax, ...
    shockssequence, irfshock, nperiods, maxiter);

Unpacking and plotting results

After the solver returns, the script unpacks the output arrays into named workspace variables using a loop over the endogenous variable list, then calls makechart to plot piecewise-linear and linear paths side by side.

% Unpack all endogenous variables from the output matrices
for i = 1:Mbase_.endo_nbr
  eval([Mbase_.endo_names{i}, '_l = zdatalinear(:,i);']);
  eval([Mbase_.endo_names{i}, '_p = zdatapiecewise(:,i);']);
  eval([Mbase_.endo_names{i}, '_ss = zdatass(i);']);
end

% Build chart inputs
titlelist  = char('c (consumption)','k (capital)','i (investment)', ...
                   'a (tfp)','lambdak (multiplier)');
ylabels    = char('Percent','Percent','Percent','Percent','Level');
legendlist = cellstr(char('Piecewise Linear','Linear'));

line1 = 100*[c_p, k_p, i_p, a_p, (lambdak_p + lambdak_ss)/100];
line2 = 100*[c_l, k_l, i_l, a_l, (lambdak_l + lambdak_ss)/100];

makechart(titlelist, legendlist, '', ylabels, line1, line2)