Occbin:Dynare identifies the wrong equilibrium condition

Hello everyone, I’ve recently encountered some issues while trying to implement occbin. In my model, there is an equilibrium condition:

E_t \leq \bar{E}

According to my model assumptions, when E_t < \bar{E}, its price \tau should be 0; when E_t = \bar{E}, its price should be 0.097890 at the steady state according to my calculations.

I wanted the system to initially be in a binding steady state, so I input the corresponding steady state into the steady_state_model block. However, it seems that Dynare considers my constraint to be slack, thus using the corresponding equilibrium condition, which leads to non-zero residuals for the equations.

I really don’t know what went wrong and hope someone can help me. Thanks in advance!

Here are my complete code and mod file, as well as the equation residuals:

// Define Endogenous Variables
var y           ${Y}$               (long_name='Output')
    c           ${C}$               (long_name='Consumption')
    i           ${I}$               (long_name='Investment')
    k           ${K}$               (long_name='Capital Stock')
    k_star      ${K^*}$             (long_name='Capital Decision')
    n           ${N}$               (long_name='Labor')
    w           ${W}$               (long_name='Wage')
    rk          ${R^K}$             (long_name='Capital Rental Rate')
    r           ${R}$               (long_name='Interest Rate')
    m           ${M}$               (long_name='Pollution Stock')
    e           ${E}$               (long_name='Emission')
    u           ${u}$               (long_name='Abatement Effort')
    z           ${Z}$               (long_name='Abatement Cost')
    tau         ${P^E}$             (long_name='Permit Price')
    a           ${A}$               (long_name='Technology')
    Delta       ${\Delta}$          (long_name='Disaster Size')
    x           ${x}$               (long_name='Disaster Binary Variable')
    lny         ${\ln(y)}$          (long_name='Log Form Output')
    lnc         ${\ln(c)}$          (long_name='Log Form Consumption')
    lni         ${\ln(i)}$          (long_name='Log Form Investment')
    lnk         ${\ln(k)}$          (long_name='Log Form Capital')
    lnks        ${\ln(k^*)}$        (long_name='Log Form Capital Decision')
;

//define predetermined variable
predetermined_variables x;

// Define Exogenous Variables
varexo  ex       ${\varepsilon^x}$               (long_name='Disaster Shock')
        ed       ${\varepsilon^{\Delta}$         (long_name='Disaster Risk Shock')
;

// Define Parameters
parameters 
    beta        ${\beta}$           (long_name='Subjective Discount Factor')
    alpha       ${\alpha}$          (long_name='Capital Share')
    phi         ${\varphi}$         (long_name='Inverse Frisch Elasticity')
    delta       ${\delta}$          (long_name='Capital Depreciation Rate')
    deltam      ${\delta_M}$        (long_name='Pollution Decay Rate')
    theta       ${\theta}$          (long_name='Disaster Probability')   
    Delta_ss    ${\Delta}$          (long_name='Steady State Value of Disaster Size')
    theta1      ${\theta_1}$        (long_name='Abatement Cost Function Efficient')
    theta2      ${\theta_2}$        (long_name='Abatement Cost Function Efficient') 
    gamma1      ${\gamma_1}$        (long_name='Emission Function Efficient')  
    d0          ${d_0}$             (long_name='Pollution Externality Coefficient')
    d1          ${d_1}$             (long_name='Pollution Externality Coefficient')
    d2          ${d_2}$             (long_name='Pollution Externality Coefficient')
    ga          ${g_A}$             (long_name='Growth Rate')    
    rho         ${\rho}$            (long_name='Disaster Risk Smoothing Parameter')
    erow        ${E^{row}}$         (long_name='Emission rof World')
    chi         ${\chi}$            (long_name='Disutility Scale')
    ebar        ${\bar{E}}$         (long_name='Permit Quota')
;

//Calibration
beta        = 0.99;        // Subjective Discount Factor
alpha       = 0.3;        // Capital Share
phi         = 3;         // Inverse Frisch Elasticity
delta       = 0.025;        // Capital Depreciation Rate
deltam      = 0.0021;        // Pollution Decay Rate
Delta_ss    = 0.43;         // Disaster Size
theta       = 0.009;         // Disaster Probability
gamma1      = 0.601;         // Emission Function Efficiency
theta1      = 0.18;         // Abatement Cost Function Efficiency
theta2      = 2.8;         // Abatement Cost Function Efficiency
d0          = 1.4e-3;        // Parameter of the damage function
d1          = -6.672e-6;     // Parameter of the damage function
d2          = 1.465e-8;      // Parameter of the damage function
rho         = 0.8; 
ga          = 0.005;
chi         = 14.4300083796; // Disutility Scale 
erow        = 1.0653773496;   // emission rest of the world 
ebar        = 0.70421477 * 0.7; 

model;  
[name='Disaster Size ']
    log(Delta) = (1-rho)*log(Delta_ss) + rho*log(Delta(-1)) + ed;

[name='Disaster']
    x(+1) = theta + ex;

[name='Growth Rate']
    a = exp(ga + x * log(1-Delta));

[name='Output Equation']
    y = k^alpha * (n)^(1-alpha) * (1 - (d0 + d1 * m + d2 * m^2));

[name='Capital Disaster']
    k * a = k_star(-1) * exp(x * log(1-Delta));

[name='Market Clearing']
    y = c + i + z;

[name='Capital Accumulation']
    k_star = (1-delta) * k + i;

[name='Optimal Capital']
    1 = beta * (c / c(+1)) * exp(x(+1) * log(1-Delta)) * (1-delta + rk(+1));
  
[name='Consumption Euler']
    1 = beta * (c / c(+1)) * r;

[name='Labor Supply']
    chi * n^phi * c = w;

[name='Emission Function']
    e = (1 - u) * gamma1 * y;

[name='Abatement Cost']
    z = y * theta1 * u^theta2;

[name='Pollution Stock Dynamics']
    m = (1 - deltam) * m(-1) / a + e + erow;

[name='Labor Demand']
    w = (1 - alpha) * (y / n);

[name='Capital Demand']
    rk = alpha * (y / k);

[name='Abatement Effort']
    u = (tau * gamma1 / (theta1 * theta2))^(1 / (theta2-1));

[name='Log Form']
    lny = log(y);
    lnc = log(c);
    lnk = log(k);
    lni = log(i);
    lnks = log(k_star);

[name='Emission Dynamics', bind='Cap']
    e = ebar;

[name='Emission Dynamics', relax='Cap']
    tau = 0;

end;

steady_state_model;
    y = 1.1771653161;
    c = 0.9061690278;
    i = 0.2634963021;
    k = 8.7795458217;
    k_star = 8.8235534782;
    n = 0.5010313778;
    w = 1.6446389544;
    rk = 0.0402241303;
    r = 1.0101010101;
    a = 0.9999409315;
    m = 763.4923833657;
    e = ebar;
    u = 0.3032271177;
    z = 0.0074999863;
    tau = 0.0978900460;
    x = theta;
    Delta = Delta_ss;
    lny = log(y);
    lnc = log(c);
    lnk = log(k);
    lnks = log(k_star);
    lni = log(i);
end;

write_latex_dynamic_model;
write_latex_static_model;

resid;
steady;

occbin_constraints;
    name 'Cap'; 
    bind e >= ebar;
    relax e < ebar;
end;

occbin_setup;


 Residuals of the static equations:
Equation number  1: Disaster Size            :       0.000000
Equation number  2: Disaster                 :       0.000000
Equation number  3: Growth Rate              :       0.000000
Equation number  4: Output Equation          :       0.000000
Equation number  5: Capital Disaster         :       0.000000
Equation number  6: Market Clearing          :       0.000000
Equation number  7: Capital Accumulation     :       0.000000
Equation number  8: Optimal Capital          :       0.000000
Equation number  9: Consumption Euler        :       0.000000
Equation number 10: Labor Supply             :       0.000000
Equation number 11: Emission Dynamics        :       0.097890
Equation number 12: Emission Function        :       0.000000
Equation number 13: Abatement Cost           :       0.000000
Equation number 14: Pollution Stock Dynamics :       0.000000
Equation number 15: Labor Demand             :       0.000000
Equation number 16: Capital Demand           :       0.000000
Equation number 17: Abatement Effort         :       0.000000
Equation number 18: Log Form                 :       0.000000
Equation number 19: lnc                      :       0.000000
Equation number 20: lnk                      :       0.000000
Equation number 21: lni                      :       0.000000
Equation number 22: lnks                     :       0.000000

What is the baseline regime, i.e. the regime where the system settles without any shocks in the long-run? The problem you describe sounds like a case where the non-binding constraint is actually the baseline.

Thank you, Professor, for your response!

I will describe my problem in more detail. My model is a DSGE model with environmental factors (E-DSGE). E_t represents pollution, which is positively correlated with output. Generally, the government implements an environmental policy to limit total pollution emissions, i.e.,
E_t \leq \bar{E}.

In the past, most approaches assumed that this constraint is always binding, i.e., E_t = \bar{E} , and solved the model under this assumption.

However, theoretically, when the economy experiences a significantly negative TFP shock, firms may substantially reduce output, leading to a considerable decrease in the demand for pollution emissions, and thus the pollution constraint might be slack, i.e., E_t < \bar{E} . According to economic theory, the price of pollution should be zero in this case.

Therefore, I want to simulate a scenario where the economy initially is at the steady state corresponding to E_t = \bar{E} . Since the constraint is binding at this point, the corresponding price of pollution is \tau > 0. Then, the economy experiences a significantly negative TFP shock, and I want to allow E_t < \bar{E} to occur during the recovery process to the initial steady state.

My code is as follows:

occbin_constraints;

    name 'Cap';   
    
    bind e >= ebar;
    
    relax e < ebar;
    
end; 
[name='Emission Dynamics', bind='Cap']
    e = ebar; 

[name='Emission Dynamics', relax='Cap']
    tau = 0;

So my baseline regime should correspond to the state E_t = \bar{E} .

Supplement:
When I changed the state corresponding to the equilibrium condition, I did obtain a steady-state solution, but the model did not converge.

[name='Emission Dynamics', relax='Cap']
    e = ebar; 

[name='Emission Dynamics', bind='Cap']
    tau = 0;

Occbin: Simulation did not converge, increase maxit or check_ahead_periods.

出错 occbin.solver (第 91 行)
    print_info(error_flag, options_.noprint, options_)

出错 simul_cap.driver (第 703 行)
[oo_.dr, oo_.occbin.simul]= occbin.solver(M_, options_, oo_.dr , oo_.steady_state, oo_.exo_steady_state, oo_.exo_det_steady_state);

出错 dynare (第 310 行)
    evalin('base',[fname '.driver']);

Firstly, I do not understand why this method yields a steady state, which contradicts my logic.

Secondly, my model runs normally without considering occbin. Once occbin is considered, the model fails to converge. How should I resolve this?

From what I understand, the relax regime is actually the one where the constraint is binding, because the occasionally binding constraint is the temporary relaxation. The long-run steady state is therefore the one where the constraint is binding. tests/occbin/model_borrcon/borrcon.mod · master · Dynare / dynare · GitLab is an example of this.
Can you provide the code for that case, which is not converging?

Dear Professor,
I found the mistake in my model that might have caused the non-convergence issue. But now I have a new question. After reviewing the example you shared, I believe I might have misunderstood the syntax of occbin. In your example, the borrowing constraint is:

B_t \leq MY_t

Assuming the multiplier corresponding to this constraint is \lambda_t, according to my understanding, when \lambda_t > 0, the borrowing constraint should be binding, i.e., B_t = MY_t based on the complementary slackness condition. When B_t < MY_t, \lambda_t = 0. Therefore, the corresponding code I think should be:

occbin_constraints;
name ‘borrcon’;
bind lb>0;
relax b<M*y;
end;
[name = ‘borrowing’, bind=‘borrcon’]
b = M*y;
[name = ‘borrowing’, relax=‘borrcon’]
lb = 0;

Unfortunately, this is exactly the opposite of your code. Where is my misunderstanding? I hope you can point it out.

The example above is one where in the baseline regime the constraint is always binding. Thus, the setup is kind of reversed. The occasionally binding constraint for implementation purposes becomes the borrowing constraint temporarily not binding, i.e. the multiplier being 0. This example shows that the OccBin implementation is about appropriately defining the two regimes, even if the “occasionally binding” part may be reversed compared to the economic intuition. I pointed you to this example because your problem is similar with the constraint binding in the baseline regime.