Overhaul Occbin interface

A two constraint-example should look like:

var a, c, i, k, lambdak;    
varexo erra;
parameters ALPHA, DELTAK, BETA, GAMMAC, RHOA, PHI, PSI;

BETA=0.96;
ALPHA=0.33;
DELTAK=0.10;
GAMMAC=2;
RHOA = 0.9;
PHI = 0.975;
PSI = 5;        % adjustment cost for capital if investment is negative

model;

# zkss = ((1/BETA-1+DELTAK)/ALPHA)^(1/(ALPHA-1));
# zcss = -DELTAK*zkss + zkss^ALPHA;
# ziss = DELTAK*zkss;
# zuss = (zcss^(1-GAMMAC)-1)/(1-GAMMAC);
# zvss = zuss/(1-BETA);

/////////////////////////////////////////////////////////////////
// 1.
[name='Euler', bind = 'INEG'] 
-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)*lambdak(1);

[name='Euler', relax = 'INEG']
-exp(c)^(-GAMMAC) + BETA*exp(c(1))^(-GAMMAC)*(1-DELTAK)+ALPHA*exp(a(1))*exp(k)^(ALPHA-1))= 
  -lambdak+BETA*(1-DELTAK)*lambdak(1);

// 2.
[name='Budget constraint',bind = 'INEG'] 
exp(c)+exp(k)-(1-DELTAK)*exp(k(-1))+PSI*(exp(k)/exp(k(-1))-1)^2=exp(a)*exp(k(-1))^(ALPHA);
[name='Budget constraint',relax = 'INEG']
exp(c)+exp(k)-(1-DELTAK)*exp(k(-1))=exp(a)*exp(k(-1))^(ALPHA);

// 3.
[name='LOM capital']
exp(i) = exp(k)-(1-DELTAK)*exp(k(-1));

// 4.
[name='investment',bind='IRR,INEG']
(i - log(PHI*ziss)) = 0;
[name='investment',relax='IRR']
lambdak=0;
[name='investment',bind='IRR',relax='INEG']
(i - log(PHI*ziss)) = 0;

// 5. 
[name='LOM TFP']
a = RHOA*a(-1)+erra;


end;

occbin_constraints;
name 'IRR';
bind i<PHI-1;
relax lambdak<0;
error_bind EXPRESSION;
error_relax EXPRESSION;
name 'INEG';
bind i<-0.000001;
end;

steady_state_model;
kss = ((1/BETA-1+DELTAK)/ALPHA)^(1/(ALPHA-1));
css = -DELTAK*kss +kss^ALPHA;
iss = DELTAK*kss;
k = log(kss);
c = log(css);
i = log(iss);
lambdak = 0;
a=0;
end;


shocks;
  var erra; stderr 0.015;
end;

steady;

  
shocks(surprise,overwrite);
var erra;
periods 1:9, 10, 50, 90, 130, 131:169;
values -0.0001, -0.01,-0.02, 0.01, 0.02, 0;
end;

occbin_setup;

occbin_solver(simul_periods=200,simul_maxit=200,simul_curb_retrench,simul_check_ahead_periods=200);
occbin_write_regimes(filename='test',periods=[1:100]);

This involves:

• Add a new block shocks(surprise)-block for setting the values of surprise shocks appearing in each period.

  shocks(surprise,overwrite);
  var NAME;
  periods INTEGER[:INTEGER] [[,] INTEGER[:INTEGER]]...;
  values DOUBLE | (EXPRESSION) [[,] DOUBLE | (EXPRESSION) ]...;
  end;

  that should map into M_.surprise_shocks.
• A new occbin_constraints-block of the form

  occbin_constraints;
  name 'STRING'; bind EXPRESSION; [relax EXPRESSION;] [error_bind EXPRESSION;] [error_relax EXPRESSION;]
  end;

  that fills the structure M_.occbin.constraint. The relax-part specifying the condition when the solution returns to the reference regime is optional. Generally, it will be simply the complement to bind.
• Remove the existing occbin-model tag and instead trigger occbin.initialize if occbin_constraints is present.
• Parse the contents written to M_.occbin.constraint (the job currently done in occbin.process_constraint and occbin.process_error_constraint (called by occbin.initialize)
• Based on the parsed content, create the +dname\eval_difference.m file currently created in occbin.initialize.
• Adjust the constraints to refer to the nonlinear instead of the linear and demeaned model by adding the steady state to all endogenous variables in the occbin_difference.m-file. Instead of

  occbin_constraints;
  name 'relax_borrcon'; bind lb<-STEADY_STATE(lb);
  end;

  to get bind lb<0 we would use

  occbin_constraints;
  name 'relax_borrcon'; bind lb<0;
  end;

  For this to work, we need to replace the current

  binding.constraint_1 = zdatalinear(:,4)<(-(steady_state(4)));

  internally by

  binding.constraint_1 = zdatalinear(:,4)+steady_state(4)<0;

• Construct the dynamic and static-files based on the present equation tags.
  • Require name-tags when bind/relax-tags are present
  • Base the static-file on the reference regime, i.e. where all constraints are relax
  • Introduce an internal binding constraint-indicating parameter based on the name of the constraints, e.g. IRR_bind, initialized at 0 and keeping track of its index in M_.occbin.constraint(k).pswitch_index
  • Construct the dynamic model to contain equations for each name-tagged equation that are a sum of the occbin-tagged ones, premultiplied with the switch-indicator. For example

  [name='investment',bind='IRR,INEG']
  (i - log(PHI*ziss)) = 0;
  [name='investment',relax='IRR']
  lambdak=0;
  [name='investment',bind='IRR',relax='INEG']
  (i - log(PHI*ziss)) = 0;

  should results in

  IRR_bind*INEG_bind*(i - log(PHI*ziss))+(1-IRR_bind)*(lambdak)+IRR_bind*(1-INEG_bind)*(i - log(PHI*ziss))