Occbin simple model consistent ss

Hi there,

I am trying to run the occbin tool on a simple model with occasionally binding constraint on one of the factor inputs, but I am not sure how to set the steady state that’s consistent with the model. I am starting with this very simple set-up, does anyone have any feedback on what’s going on?

//===============================================
// Title: One-Sector Model with OccBin Constraint
//        on Materials M (Binding and Non-Binding Versions)
// Y = AK^alpha M^(1-alpha); M <= bar{M}
//===============================================

//-----------------------------------------------
// 1. VARIABLES AND EXOGENOUS SHOCK
//-----------------------------------------------
var 
    C      (long_name='Consumption')
    K      (long_name='Capital Stock')
    Y      (long_name='Output')
    I      (long_name='Investment')
    M      (long_name='Materials Used')
    A      (long_name='Productivity')
    lambdak (long_name='Lagrangian on constraint');

varexo erra;  // Productivity shock

//-----------------------------------------------
// 2. PARAMETERS
//-----------------------------------------------
parameters beta alpha delta rho sigma M_bar;

beta   = 0.96;       // Discount factor
alpha  = 0.33;       // Capital share
delta  = 0.025;      // Depreciation rate
rho    = 0.95;       // AR(1) parameter for productivity
sigma  = 0.01;       // Std dev of productivity shock
M_bar  = 1.5;        // Upper bound on Materials

//-----------------------------------------------
// 3. MODEL BLOCK (Binding Case)
//-----------------------------------------------
model;
    [name='Regime1_Binding']
    // (1) Cobb-Douglas production function
    Y = A * K^alpha * M^(1 - alpha);

    // (2) Capital accumulation equation
    K = (1 - delta) * K(-1) + I;

    // (3) Goods market clearing condition
    Y = C + I;

    // (4) Euler equation
    1/C = beta * (1/C(+1)) * ( alpha * Y(+1)/K(+1) + (1 - delta) );

    // (5) Market equilibrium for materials
    (1 - alpha) * Y / M = C * lambdak;

    // (6) Productivity AR(1) process
    log(A) = rho * log(A(-1)) + erra;

     // (7) MATERIALS SUPPLY CONSTRAINT (Binding Version)
     [name='Materials constraint', bind = 'shift']
     M = M_bar;
     [name='Materials constraint', relax = 'shift']
     lambdak = 0;

end;


//-----------------------------------------------
// 4. INITIAL STEADY STATE GUESS
//-----------------------------------------------
initval;
    A = 1;                  // Productivity normalized
    M = M_bar;              // Materials at upper bound (binding case)
    erra = 0;               // No shock in steady state

    // Solve for steady-state capital using Euler equation
    K  = M*(((1/beta) - 1 + delta) / (alpha * A))^(1 / (alpha - 1)); 
    
    // Output from Cobb-Douglas production function
    Y = A * K^alpha * M^(1 - alpha);    
    
    // Investment equals depreciation in steady state
    I = delta * K;          
    
    // Consumption from goods market clearing
    C = Y - I;              
    
    // Solve for lambda_k using materials market equilibrium
    lambdak = (1 - alpha) * Y / (M * C);
end;

check;
steady;

//-----------------------------------------------
// 5. OCCBIN CONSTRAINTS DEFINITION
//-----------------------------------------------
occbin_constraints;
    name 'shift'; bind   M < M_bar;        // The constraint binds when M exceeds M_bar
end;

//-----------------------------------------------
// 6. SHOCKS
//-----------------------------------------------
shocks(surprise,overwrite);
    var erra;  % Specify variable to shock
    periods 1 2 3 4 5 6;  % Specify when to give the shock in the simulation.
    values 0.02 -0.01 0.015 -0.005 0 0;
end;

//-----------------------------------------------
// 7. OCCBIN SETUP & SIMULATION
//-----------------------------------------------
occbin_setup; 
occbin_solver;

// Compute and plot the IRFs using OccBin
occbin_graph;

The steady state needs to be provided for the relax regime. Why do you set lambdak = 0; in that regime but use

    lambdak = (1 - alpha) * Y / (M * C);

for the initial value?