diff --git a/matlab/rotated_slice_sampler.m b/matlab/rotated_slice_sampler.m
new file mode 100644
index 0000000000000000000000000000000000000000..66908d77d44264a54693061b5a7cb738bb64231a
--- /dev/null
+++ b/matlab/rotated_slice_sampler.m
@@ -0,0 +1,184 @@
+function [theta, fxsim, neval] = rotated_slice_sampler(objective_function,theta,thetaprior,sampler_options,varargin)
+% ----------------------------------------------------------
+% ROTATED SLICE SAMPLER - with stepping out (Neal, 2003)
+% extension of the orthogonal univarite sampler (slice_sampler.m)
+% copyright M. Ratto (European Commission)
+%
+% objective_function(theta,varargin): -log of any unnormalized pdf
+% with varargin (optional) a vector of auxiliaty parameters
+% to be passed to f( ).
+% ----------------------------------------------------------
+%
+% INPUTS
+% objective_function: objective function (expressed as minus the log of a density)
+% theta: last value of theta
+% thetaprior: bounds of the theta space
+% sampler_options: posterior sampler options
+% varargin: optional input arguments to objective function
+%
+% OUTPUTS
+% theta: new theta sample
+% fxsim: value of the objective function for the new sample
+% neval: number of function evaluations
+%
+% SPECIAL REQUIREMENTS
+% none
+
+% 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/>.
+
+theta=theta(:);
+npar = length(theta);
+neval = zeros(npar,1);
+W1=[];
+if isfield(sampler_options,'WR'),
+ W1 = sampler_options.WR;
+end
+if ~isempty(sampler_options.mode),
+ mm = sampler_options.mode;
+ n = length(mm);
+ for j=1:n,
+ distance(j)=sqrt(sum((theta-mm(j).m).^2));
+ end
+ [m, im] = min(distance);
+
+ r=im;
+ V1 = mm(r).m;
+ jj=0;
+ for j=1:n,
+ if j~=r,
+ jj=jj+1;
+ tmp=mm(j).m-mm(r).m;
+ %tmp=mm(j).m-theta;
+ V1(:,jj)=tmp/norm(tmp);
+ end
+ end
+ resul=randperm(n-1,n-1);
+ V1 = V1(:,resul);
+
+ %V1 = V1(:, randperm(n-1));
+% %d = chol(mm(r).invhess);
+% %V1 = transpose(feval(sampler_options.proposal_distribution, transpose(mm(r).m), d, npar));
+%
+% V1=eye(npar);
+% V1=V1(:,randperm(npar));
+% for j=1:2,
+% V1(:,j)=mm(r(j)).m-theta;
+% V1(:,j)=V1(:,j)/norm(V1(:,j));
+% end
+% % Gram-Schmidt
+% for j=2:npar,
+% for k=1:j-1,
+% V1(:,j)=V1(:,j)-V1(:,k)'*V1(:,j)*V1(:,k);
+% end
+% V1(:,j)=V1(:,j)/norm(V1(:,j));
+% end
+% for j=1:n,
+% distance(j)=sqrt(sum((theta-mm(j).m).^2));
+% end
+% [m, im] = min(distance);
+% if im==r,
+% fxsim=[];
+% return,
+% else
+% theta1=theta;
+% end
+else
+ V1 = sampler_options.V1;
+end
+npar=size(V1,2);
+
+for it=1:npar,
+ theta0 = theta;
+ neval(it) = 0;
+ xold = 0;
+ % XLB = thetaprior(3);
+ % XUB = thetaprior(4);
+ tb=sort([(thetaprior(:,1)-theta)./V1(:,it) (thetaprior(:,2)-theta)./V1(:,it)],2);
+ XLB=max(tb(:,1));
+ XUB=min(tb(:,2));
+ if isempty(W1),
+ W = (XUB-XLB); %*0.8;
+ else
+ W = W1(it);
+ end
+
+ % -------------------------------------------------------
+ % 1. DRAW Z = ln[f(X0)] - EXP(1) where EXP(1)=-ln(U(0,1))
+ % THIS DEFINES THE SLICE S={x: z < ln(f(x))}
+ % -------------------------------------------------------
+
+ fxold = -feval(objective_function,theta,varargin{:});
+ %I have to be sure that the rotation is for L,R or for Fxold, theta(it)
+ neval(it) = neval(it) + 1;
+ Z = fxold + log(rand(1,1));
+ % -------------------------------------------------------------
+ % 2. FIND I=(L,R) AROUND X0 THAT CONTAINS S AS MUCH AS POSSIBLE
+ % STEPPING-OUT PROCEDURE
+ % -------------------------------------------------------------
+ u = rand(1,1);
+ L = max(XLB,xold-W*u);
+ R = min(XUB,L+W);
+
+ %[L R]=slice_rotation(L, R, alpha);
+ while(L > XLB)
+ xsim = L;
+ theta = theta0+xsim*V1(:,it);
+ fxl = -feval(objective_function,theta,varargin{:});
+ neval(it) = neval(it) + 1;
+ if (fxl <= Z)
+ break;
+ end
+ L = max(XLB,L-W);
+ end
+ while(R < XUB)
+ xsim = R;
+ theta = theta0+xsim*V1(:,it);
+ fxr = -feval(objective_function,theta,varargin{:});
+ neval(it) = neval(it) + 1;
+ if (fxr <= Z)
+ break;
+ end
+ R = min(XUB,R+W);
+ end
+ % ------------------------------------------------------
+ % 3. SAMPLING FROM THE SET A = (I INTERSECT S) = (LA,RA)
+ % ------------------------------------------------------
+ fxsim = Z-1;
+ while (fxsim < Z)
+ u = rand(1,1);
+ xsim = L + u*(R - L);
+ theta = theta0+xsim*V1(:,it);
+ fxsim = -feval(objective_function,theta,varargin{:});
+ neval(it) = neval(it) + 1;
+ if (xsim > xold)
+ R = xsim;
+ else
+ L = xsim;
+ end
+ end
+end
+
+% if ~isempty(sampler_options.mode),
+% dist1=sqrt(sum((theta-mm(r).m).^2));
+% if dist1>distance(r),
+% theta=theta1;
+% fxsim=[];
+% end
+% end
+end
+
diff --git a/matlab/slice_sampler.m b/matlab/slice_sampler.m
new file mode 100644
index 0000000000000000000000000000000000000000..e454d59365da8030f5b871d5d8ee780c9e7628b4
--- /dev/null
+++ b/matlab/slice_sampler.m
@@ -0,0 +1,123 @@
+function [theta, fxsim, neval] = slice_sampler(objective_function,theta,thetaprior,sampler_options,varargin)
+% function [theta, fxsim, neval] = slice_sampler(objective_function,theta,thetaprior,sampler_options,varargin)
+% ----------------------------------------------------------
+% UNIVARIATE SLICE SAMPLER - stepping out (Neal, 2003)
+% W: optimal value in the range (3,10)*std(x)
+% - see C.Planas and A.Rossi (2014)
+% objective_function(theta,varargin): -log of any unnormalized pdf
+% with varargin (optional) a vector of auxiliaty parameters
+% to be passed to f( ).
+% ----------------------------------------------------------
+%
+% INPUTS
+% objective_function: objective function (expressed as minus the log of a density)
+% theta: last value of theta
+% thetaprior: bounds of the theta space
+% sampler_options: posterior sampler options
+% varargin: optional input arguments to objective function
+%
+% OUTPUTS
+% theta: new theta sample
+% fxsim: value of the objective function for the new sample
+% neval: number of function evaluations
+%
+% SPECIAL REQUIREMENTS
+% none
+
+% 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/>.
+
+if sampler_options.rotated %&& ~isempty(sampler_options.V1),
+ [theta, fxsim, neval] = rotated_slice_sampler(objective_function,theta,thetaprior,sampler_options,varargin{:});
+ if isempty(sampler_options.mode), % jumping
+ return,
+ else
+ nevalR=sum(neval);
+ end
+end
+
+theta=theta(:);
+npar = length(theta);
+W1 = sampler_options.W1;
+neval = zeros(npar,1);
+
+for it=1:npar,
+ neval(it) = 0;
+ W = W1(it);
+ xold = theta(it);
+ % XLB = thetaprior(3);
+ % XUB = thetaprior(4);
+ XLB = thetaprior(it,1);
+ XUB = thetaprior(it,2);
+
+
+ % -------------------------------------------------------
+ % 1. DRAW Z = ln[f(X0)] - EXP(1) where EXP(1)=-ln(U(0,1))
+ % THIS DEFINES THE SLICE S={x: z < ln(f(x))}
+ % -------------------------------------------------------
+ fxold = -feval(objective_function,theta,varargin{:});
+ neval(it) = neval(it) + 1;
+ Z = fxold + log(rand(1,1));
+ % -------------------------------------------------------------
+ % 2. FIND I=(L,R) AROUND X0 THAT CONTAINS S AS MUCH AS POSSIBLE
+ % STEPPING-OUT PROCEDURE
+ % -------------------------------------------------------------
+ u = rand(1,1);
+ L = max(XLB,xold-W*u);
+ R = min(XUB,L+W);
+ while(L > XLB)
+ xsim = L;
+ theta(it) = xsim;
+ fxl = -feval(objective_function,theta,varargin{:});
+ neval(it) = neval(it) + 1;
+ if (fxl <= Z)
+ break;
+ end
+ L = max(XLB,L-W);
+ end
+ while(R < XUB)
+ xsim = R;
+ theta(it) = xsim;
+ fxr = -feval(objective_function,theta,varargin{:});
+ neval(it) = neval(it) + 1;
+ if (fxr <= Z)
+ break;
+ end
+ R = min(XUB,R+W);
+ end
+ % ------------------------------------------------------
+ % 3. SAMPLING FROM THE SET A = (I INTERSECT S) = (LA,RA)
+ % ------------------------------------------------------
+ fxsim = Z-1;
+ while (fxsim < Z)
+ u = rand(1,1);
+ xsim = L + u*(R - L);
+ theta(it) = xsim;
+ fxsim = -feval(objective_function,theta,varargin{:});
+ neval(it) = neval(it) + 1;
+ if (xsim > xold)
+ R = xsim;
+ else
+ L = xsim;
+ end
+ end
+
+end
+
+if sampler_options.rotated && ~isempty(sampler_options.mode), % jumping
+ neval=sum(neval)+nevalR;
+end