diff --git a/matlab/occbin/evaluate_model.m b/matlab/occbin/evaluate_model.m
new file mode 100644
index 0000000000000000000000000000000000000000..f01d58a08df113be9ee59b9a1311ca612e7ed458
--- /dev/null
+++ b/matlab/occbin/evaluate_model.m
@@ -0,0 +1,19 @@
+function [r,g1,g2,g3] = evaluate_model(z,x,M,ss)
+
+ll = M.lead_lag_incidence';
+y = z(find(ll(:)));
+
+switch nargout
+  case 1
+    r = feval([M.fname '_dynamic'],y,x, ...
+              M.params, ss, 1);
+  case 2
+    [r,g1] = feval([M.fname '_dynamic'],y,x, ...
+                    M.params, ss, 1);
+  case 3
+    [r,g1,g2] = feval([M.fname '_dynamic'],y,x, ...
+                    M.params, ss, 1);
+  case 4
+    [r,g1,g2,g3] = feval([M.fname '_dynamic'],y,x, ...
+                         M.params, ss, 1);
+end
\ No newline at end of file
diff --git a/matlab/occbin/get_coef.m b/matlab/occbin/get_coef.m
new file mode 100644
index 0000000000000000000000000000000000000000..f2dbde5a619dce0a95c22c304f09b4ab80258019
--- /dev/null
+++ b/matlab/occbin/get_coef.m
@@ -0,0 +1,26 @@
+function [coef_y,coef_u] = get_coef(jacobian,M)
+
+ll = M.lead_lag_incidence;
+endo_nbr = M.endo_nbr;
+coef_y = zeros(endo_nbr,3*endo_nbr);
+coef_u = zeros(endo_nbr,M.exo_nbr);
+
+if M.maximum_lag > 0
+    [junk,c1,c2] = find(ll(1,:));
+    coef_y(:,c1) = jacobian(:,c2);
+    [junk,c1,c2] = find(ll(2,:));
+    coef_y(:,c1+endo_nbr) = jacobian(:,c2);
+    if M.maximum_lead > 0
+            [junk,c1,c2] = find(ll(3,:));
+            coef_y(:,c1+2*endo_nbr) = jacobian(:,c2);
+    end
+else
+    [junk,c1,c2] = find(ll(1,:));
+    coef_y(:,c1+endo_nbr) = jacobian(:,c2);
+    if M.maximum_lead > 0
+            [junk,c1,c2] = find(ll(2,:));
+            coef_y(:,c1+2*endo_nbr) = jacobian(:,c2);
+    end
+end
+
+coef_u = jacobian(:,max(ll(end,:))+1:end);
\ No newline at end of file
diff --git a/matlab/occbin/get_complementarity_conditions.m b/matlab/occbin/get_complementarity_conditions.m
new file mode 100644
index 0000000000000000000000000000000000000000..cf5a40866c68ec71ca77295ed83dd85b524949b9
--- /dev/null
+++ b/matlab/occbin/get_complementarity_conditions.m
@@ -0,0 +1,72 @@
+function [lb,ub,eq_index] = get_complementarity_conditions(M,ramsey_policy)
+
+% Copyright (C) 2014 Dynare Team
+%
+% This file is part of Dynare.
+%
+% Dynare is free software: you can redistribute it and/or modify
+% it under the terms of the GNU General Public License as published by
+% the Free Software Foundation, either version 3 of the License, or
+% (at your option) any later version.
+%
+% Dynare is distributed in the hope that it will be useful,
+% but WITHOUT ANY WARRANTY; without even the implied warranty of
+% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+% GNU General Public License for more details.
+%
+% You should have received a copy of the GNU General Public License
+% along with Dynare.  If not, see <http://www.gnu.org/licenses/>.
+
+ub = inf(M.endo_nbr,1);
+lb = -ub;
+eq_index = (1:M.endo_nbr)';
+if ramsey_policy
+    if isfield(M,'ramsey_model_constraints')
+        rc = M.ramsey_model_constraints;
+        for i = 1:length(rc)
+            switch rc{i}{2}
+              case {'>','>='}
+                lb(rc{i}{1}) = eval(rc{i}{3});
+              case {'<','<='}
+                ub(rc{i}{1}) = eval(rc{i}{3});
+              otherwise
+                error('Wrong operator in get_complementarity_conditions')
+            end
+            eq_index(i) = 1;
+        end
+    end
+end
+
+etags = M.equations_tags;
+for i=1:size(etags,1)
+    if strcmp(etags{i,2},'mcp')
+        str = etags{i,3};
+        kop = strfind(etags{i,3},'<');
+        if ~isempty(kop)
+                k = find(strcmp(strtrim(str(1:kop-1)),cellstr(M.endo_names)));
+                if isempty(k)
+                    error(sprintf(['Complementarity condition %s: variable %s is ' ...
+                                   'not recognized',etags{i,3},b{1}]))
+                end
+                ub(k) = str2num(str(kop+1:end));
+                eq_index(etags{i,1}) = k;
+                eq_index(k) = etags{i,1};
+        else
+            kop = strfind(etags{i,3},'>');
+            if ~isempty(kop)
+                k = find(strcmp(strtrim(str(1:kop-1)),cellstr(M.endo_names)));
+                if isempty(k)
+                    error(sprintf(['Complementarity condition %s: variable %s is ' ...
+                                   'not recognized',etags{i},b{1}]))
+                end
+                lb(k) = str2num(str(kop+1:end));
+                eq_index(etags{i,1}) = k;
+                eq_index(k) = etags{i,1};
+            else
+                error(sprintf(['Complementarity condition %s can''t be ' ...
+                               'parsed'],etags{i,3}))
+            end
+        end
+    end
+end
+
diff --git a/matlab/occbin/get_occbin_complementarity_conditions.m b/matlab/occbin/get_occbin_complementarity_conditions.m
new file mode 100644
index 0000000000000000000000000000000000000000..29c8655872393b95058ad4830dfb2ef302651bd9
--- /dev/null
+++ b/matlab/occbin/get_occbin_complementarity_conditions.m
@@ -0,0 +1,75 @@
+function [ivar,ieq,lb,ub] = get_occbin_complementarity_conditions(M,ramsey_policy)
+
+% Copyright (C) 2015 Dynare Team
+%
+% This file is part of Dynare.
+%
+% Dynare is free software: you can redistribute it and/or modify
+% it under the terms of the GNU General Public License as published by
+% the Free Software Foundation, either version 3 of the License, or
+% (at your option) any later version.
+%
+% Dynare is distributed in the hope that it will be useful,
+% but WITHOUT ANY WARRANTY; without even the implied warranty of
+% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+% GNU General Public License for more details.
+%
+% You should have received a copy of the GNU General Public License
+% along with Dynare.  If not, see <http://www.gnu.org/licenses/>.
+
+nrow = 1;
+if ramsey_policy
+    if isfield(M,'ramsey_model_constraints')
+        rc = M.ramsey_model_constraints;
+        for i = 1:length(rc)
+            switch rc{i}{2}
+              case {'>','>='}
+                ivar(nrow) = rc{i}{1};
+                ieq(nrow) = rc{i}{1};
+                lb(nrow) = eval(rc{i}{3});
+              case {'<','<='}
+                ivar(nrow) = rc{i}{1};
+                ieq(nrow) = rc{i}{1};
+                ub(nrow) = eval(rc{i}{3});
+              otherwise
+                error('Wrong operator in get_complementarity_conditions')
+            end
+            nrow = nrow + 1;
+        end
+    end
+end
+
+etags = M.equations_tags;
+for i=1:size(etags,1)
+    if strcmp(etags{i,2},'mcp')
+        str = etags{i,3};
+        kop = strfind(etags{i,3},'<');
+        if ~isempty(kop)
+                k = find(strcmp(strtrim(str(1:kop-1)),cellstr(M.endo_names)));
+                if isempty(k)
+                    error(sprintf(['Complementarity condition %s: variable %s is ' ...
+                                   'not recognized',etags{i,3},b{1}]))
+                end
+                ivar(nrow) = k;
+                ieq(nrow) = etags{i,1};
+                ub(nrow) = eval(str(kop+1:end));
+        else
+            kop = strfind(etags{i,3},'>');
+            if ~isempty(kop)
+                k = find(strcmp(strtrim(str(1:kop-1)),cellstr(M.endo_names)));
+                if isempty(k)
+                    error(sprintf(['Complementarity condition %s: variable %s is ' ...
+                                   'not recognized',etags{i},b{1}]))
+                end
+                ivar(nrow) = k;
+                ieq(nrow) = etags{i,1};
+                lb(k) = eval(str(kop+1:end));
+            else
+                error(sprintf(['Complementarity condition %s can''t be ' ...
+                               'parsed'],etags{i,3}))
+            end
+        end
+        nrow = nrow + 1;
+    end
+end
+
diff --git a/matlab/occbin/get_occbin_constraints.m b/matlab/occbin/get_occbin_constraints.m
new file mode 100644
index 0000000000000000000000000000000000000000..6723f45423b13e99879900f0ec61413a9ccb9387
--- /dev/null
+++ b/matlab/occbin/get_occbin_constraints.m
@@ -0,0 +1,88 @@
+function [i_base,i_alt,c_base,c_alt] = get_occbin_constraints(M,steady_state,ramsey_policy)
+
+% Copyright (C) 2015 Dynare Team
+%
+% This file is part of Dynare.
+%
+% Dynare is free software: you can redistribute it and/or modify
+% it under the terms of the GNU General Public License as published by
+% the Free Software Foundation, either version 3 of the License, or
+% (at your option) any later version.
+%
+% Dynare is distributed in the hope that it will be useful,
+% but WITHOUT ANY WARRANTY; without even the implied warranty of
+% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+% GNU General Public License for more details.
+%
+% You should have received a copy of the GNU General Public License
+% along with Dynare.  If not, see <http://www.gnu.org/licenses/>.
+
+nrow = 1;
+if ramsey_policy
+    if isfield(M,'ramsey_model_constraints')
+        rc = M.ramsey_model_constraints;
+        for i = 1:length(rc)
+            switch rc{i}{2}
+              case {'>','>='}
+                ivar(nrow) = rc{i}{1};
+                ieq(nrow) = rc{i}{1};
+                lb(nrow) = eval(rc{i}{3});
+              case {'<','<='}
+                ivar(nrow) = rc{i}{1};
+                ieq(nrow) = rc{i}{1};
+                ub(nrow) = eval(rc{i}{3});
+              otherwise
+                error('Wrong operator in get_complementarity_conditions')
+            end
+            nrow = nrow + 1;
+        end
+    end
+end
+
+i_base = {};
+i_alt = {};
+etags = M.equations_tags;
+m = 1;
+base = true;
+for i=1:size(etags,1)
+    [iv,boundary,operator] = parse_constraint(etags{i,3},M.endo_names,M.params,M.param_names);
+    if strcmp(etags{i,2},'OCCBIN')
+        if base
+            i_alt{m} = 1:M.eq_nbr;
+            i_alt{m}(etags{i,1}) = [];
+            c_base{m,1} = etags{i,1};
+            c_base{m,2} = iv;
+            c_base{m,3} = boundary - steady_state(iv);
+            c_base{m,4} = operator;
+            base = false;
+        else
+            i_base{m} = 1:M.eq_nbr;
+            i_base{m}(etags{i,1}) = [];
+            c_alt{m,1} = etags{i,1};
+            c_alt{m,2} = iv;
+            c_alt{m,3} = boundary - steady_state(iv);
+            c_alt{m,4} = operator;
+            base = true;
+            m = m + 1;
+        end
+    end
+end
+if ~base
+    error('OCCBIN: constraints must come by pair')
+end
+
+function [iv,boundary,operator] = parse_constraint(str,endo_names,params,param_names)
+delim = {'<=','>=','<','>'};
+[c,operator] = strsplit(str,delim);
+operator = operator{1};
+iv = strmatch(strtrim(c{1}),endo_names);
+% try for a number
+boundary = str2num(strtrim(c{2}));
+% if not a number try for a parameter name
+if isempty(boundary)
+    k = strmatch(strtrim(c{2}),param_names);
+    if isempty(k)
+        error(['OCCBIN: illegal constraint ' str]);
+    end
+    boundary = params(k);
+end
diff --git a/matlab/occbin/get_residuals.m b/matlab/occbin/get_residuals.m
new file mode 100644
index 0000000000000000000000000000000000000000..f86d05cfb2e6802be178ff5250dfc788943a6473
--- /dev/null
+++ b/matlab/occbin/get_residuals.m
@@ -0,0 +1,10 @@
+function r = get_residuals(ivar,lb,ub,M,oo)
+
+ss = oo.steady_state;
+for i = 1:length(ivar)
+    % only one is different from zero
+    ss(ivar(i)) = lb(i) + ub(i);
+end
+oo.steady_state = ss;
+
+r = evaluate_model(M,oo);
\ No newline at end of file
diff --git a/matlab/occbin/occbin.m b/matlab/occbin/occbin.m
new file mode 100644
index 0000000000000000000000000000000000000000..ea0780d2111016114fc3da7203bde998798ff6b8
--- /dev/null
+++ b/matlab/occbin/occbin.m
@@ -0,0 +1,19 @@
+function [endo_simul,endo_simul_no_constraint,status] = occbin(M,oo,options)
+% function oo=occbin(M,oo,options) solves linear models with occasionally
+% binding constraints using OCCBIN by L. Guerrieri
+
+status = 1;
+constraint_nbr = sum(strcmp(upper(M.equations_tags(:,2)),'OCCBIN'))/2;
+
+switch(constraint_nbr)
+  case 1
+    [zdatalinear_ zdatapiecewise_ zdatass_ oobase_ ] = ...
+    solve_one_constraint(M,oo,options);
+  case 2
+    [zdatalinear_ zdatapiecewise_ zdatass_ oobase_ ] = ...
+    solve_two_constraints(M,oo,options);
+  otherwise
+    error('OCCBIN can only handle two constraints in a model')
+end
+endo_simul = zdatapiecewise_';
+endo_simul_no_constraint = zdatalinear_';
\ No newline at end of file
diff --git a/matlab/occbin/test_constraint.m b/matlab/occbin/test_constraint.m
new file mode 100644
index 0000000000000000000000000000000000000000..96e41839a2e1ef8e372e4d63f022d17a279d91ea
--- /dev/null
+++ b/matlab/occbin/test_constraint.m
@@ -0,0 +1,16 @@
+function b=test_constraint(x,constr)
+b = zeros(size(x,1),size(constr,1));
+for i=1:size(constr,1)
+    switch constr{i,4}
+      case '<'
+        b(:,i) = ~(x(:,constr{i,2}) < constr{i,3});
+      case '<='
+        b(:,i) = ~(x(:,constr{i,2}) <= constr{i,3});
+      case '>'
+        b(:,i) = ~(x(:,constr{i,2}) > constr{i,3});
+      case '>='
+        b(:,i) = ~(x(:,constr{i,2}) >= constr{i,3});
+      otherwise
+        error('OCCBIN: wrong inequality sign')
+    end
+end
\ No newline at end of file