diff --git a/matlab/AHessian.m b/matlab/AHessian.m
deleted file mode 100644
index 64225197fee159e7bf280b97d577cf9cbea9b35f..0000000000000000000000000000000000000000
--- a/matlab/AHessian.m
+++ /dev/null
@@ -1,152 +0,0 @@
-function [AHess, DLIK, LIK] = AHessian(T,R,Q,H,P,Y,DT,DYss,DOm,DH,DP,start,mf,kalman_tol,riccati_tol)
-% function [AHess, DLIK, LIK] = AHessian(T,R,Q,H,P,Y,DT,DYss,DOm,DH,DP,start,mf,kalman_tol,riccati_tol)
-%
-% computes the asymptotic hessian matrix of the log-likelihood function of
-% a state space model (notation as in kalman_filter.m in DYNARE
-% Thanks to  Nikolai Iskrev
-%
-% NOTE: the derivative matrices (DT,DR ...) are 3-dim. arrays with last
-% dimension equal to the number of structural parameters
-
-% Copyright © 2011-2017 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 <https://www.gnu.org/licenses/>.
-
-
-k = size(DT,3);                                 % number of structural parameters
-smpl = size(Y,2);                               % Sample size.
-pp   = size(Y,1);                               % Maximum number of observed variables.
-mm   = size(T,2);                               % Number of state variables.
-a    = zeros(mm,1);                             % State vector.
-Om   = R*Q*transpose(R);                        % Variance of R times the vector of structural innovations.
-t    = 0;                                       % Initialization of the time index.
-oldK = 0;
-notsteady   = 1;                                % Steady state flag.
-F_singular  = 1;
-
-lik  = zeros(smpl,1);                           % Initialization of the vector gathering the densities.
-LIK  = Inf;                                     % Default value of the log likelihood.
-if nargout > 1
-    DLIK  = zeros(k,1);                             % Initialization of the score.
-end
-AHess  = zeros(k,k);                             % Initialization of the Hessian
-Da    = zeros(mm,k);                             % State vector.
-Dv = zeros(length(mf),k);
-
-%     for ii = 1:k
-%         DOm = DR(:,:,ii)*Q*transpose(R) + R*DQ(:,:,ii)*transpose(R) + R*Q*transpose(DR(:,:,ii));
-%     end
-
-while notsteady && t<smpl
-    t  = t+1;
-    v  = Y(:,t)-a(mf);
-    F  = P(mf,mf) + H;
-    if rcond(F) < kalman_tol
-        if ~all(abs(F(:))<kalman_tol)
-            return
-        else
-            a = T*a;
-            P = T*P*transpose(T)+Om;
-        end
-    else
-        F_singular = 0;
-        iF     = inv(F);
-        K      = P(:,mf)*iF;
-        lik(t) = log(det(F))+transpose(v)*iF*v;
-
-        [DK,DF,DP1] = computeDKalman(T,DT,DOm,P,DP,DH,mf,iF,K);
-
-        for ii = 1:k
-            Dv(:,ii)   = -Da(mf,ii) - DYss(mf,ii);
-            Da(:,ii)   = DT(:,:,ii)*(a+K*v) + T*(Da(:,ii)+DK(:,:,ii)*v + K*Dv(:,ii));
-            if t>=start && nargout > 1
-                DLIK(ii,1)  = DLIK(ii,1) + trace( iF*DF(:,:,ii) ) + 2*Dv(:,ii)'*iF*v - v'*(iF*DF(:,:,ii)*iF)*v;
-            end
-        end
-        vecDPmf = reshape(DP(mf,mf,:),[],k);
-        %                     iPmf = inv(P(mf,mf));
-        if t>=start
-            AHess = AHess + Dv'*iF*Dv + .5*(vecDPmf' * kron(iF,iF) * vecDPmf);
-        end
-        a      = T*(a+K*v);
-        P      = T*(P-K*P(mf,:))*transpose(T)+Om;
-        DP     = DP1;
-    end
-    notsteady = max(max(abs(K-oldK))) > riccati_tol;
-    oldK = K;
-end
-
-if F_singular
-    error('The variance of the forecast error remains singular until the end of the sample')
-end
-
-
-if t < smpl
-    t0 = t+1;
-    while t < smpl
-        t = t+1;
-        v = Y(:,t)-a(mf);
-        for ii = 1:k
-            Dv(:,ii)   = -Da(mf,ii)-DYss(mf,ii);
-            Da(:,ii)   = DT(:,:,ii)*(a+K*v) + T*(Da(:,ii)+DK(:,:,ii)*v + K*Dv(:,ii));
-            if t>=start && nargout >1
-                DLIK(ii,1)  = DLIK(ii,1) + trace( iF*DF(:,:,ii) ) + 2*Dv(:,ii)'*iF*v - v'*(iF*DF(:,:,ii)*iF)*v;
-            end
-        end
-        if t>=start
-            AHess = AHess + Dv'*iF*Dv;
-        end
-        a = T*(a+K*v);
-        lik(t) = transpose(v)*iF*v;
-    end
-    AHess = AHess + .5*(smpl-t0+1)*(vecDPmf' * kron(iF,iF) * vecDPmf);
-    if nargout > 1
-        for ii = 1:k
-            %             DLIK(ii,1)  = DLIK(ii,1) + (smpl-t0+1)*trace( iF*DF(:,:,ii) );
-        end
-    end
-    lik(t0:smpl) = lik(t0:smpl) + log(det(F));
-    %         for ii = 1:k;
-    %             for jj = 1:ii
-    %              H(ii,jj) = trace(iPmf*(.5*DP(mf,mf,ii)*iPmf*DP(mf,mf,jj) + Dv(:,ii)*Dv(:,jj)'));
-    %             end
-    %         end
-end
-
-AHess = -AHess;
-if nargout > 1
-    DLIK = DLIK/2;
-end
-% adding log-likelihhod constants
-lik = (lik + pp*log(2*pi))/2;
-
-LIK = sum(lik(start:end)); % Minus the log-likelihood.
-                           % end of main function
-
-function [DK,DF,DP1] = computeDKalman(T,DT,DOm,P,DP,DH,mf,iF,K)
-
-k      = size(DT,3);
-tmp    = P-K*P(mf,:);
-
-for ii = 1:k
-    DF(:,:,ii)  = DP(mf,mf,ii) + DH(:,:,ii);
-    DiF(:,:,ii) = -iF*DF(:,:,ii)*iF;
-    DK(:,:,ii)  = DP(:,mf,ii)*iF + P(:,mf)*DiF(:,:,ii);
-    Dtmp        = DP(:,:,ii) - DK(:,:,ii)*P(mf,:) - K*DP(mf,:,ii);
-    DP1(:,:,ii) = DT(:,:,ii)*tmp*T' + T*Dtmp*T' + T*tmp*DT(:,:,ii)' + DOm(:,:,ii);
-end
-
-% end of computeDKalman
diff --git a/matlab/bseastr.m b/matlab/bseastr.m
deleted file mode 100644
index 1c9a888afba855338db8f0fc65c2cebc036ad7b6..0000000000000000000000000000000000000000
--- a/matlab/bseastr.m
+++ /dev/null
@@ -1,47 +0,0 @@
-function x = bseastr(s1,s2)
-
-% Copyright © 2001-2017 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 <https://www.gnu.org/licenses/>.
-
-m = size(s1,1) ;
-x = zeros(m,1) ;
-s1=upper(deblank(s1));
-s2=upper(deblank(s2));
-
-for im = 1:m
-    key = s1(im,:) ;
-    h = size(s2,1) ;
-    l = 1 ;
-    while l <= h
-        mid = round((h+l)/2) ;
-        temp = s2(mid,:) ;
-        if ~ strcmp(key,temp)
-            for i = 1:min(length(key),length(temp))
-                if temp(i) > key(i)
-                    h = mid - 1 ;
-                    break
-                else
-                    l = mid + 1 ;
-                    break
-                end
-            end
-        else
-            x(im) = mid ;
-            break
-        end
-    end
-end
diff --git a/matlab/ftest.m b/matlab/ftest.m
deleted file mode 100644
index d3a211482951e8eeaf8941ab5c68c54039f80dfc..0000000000000000000000000000000000000000
--- a/matlab/ftest.m
+++ /dev/null
@@ -1,84 +0,0 @@
-function ftest (s1,s2)
-
-% Copyright © 2001-2017 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 <https://www.gnu.org/licenses/>.
-
-global nvx nvy x y lag1
-
-if size(s1,1) ~= 2
-    error ('Spécifiez deux fichiers pour la comparaison.') ;
-end
-
-for i = 1:2
-    if ~ isempty(find(abs(s1(i,:)) == 46))
-        error ('Entrez les noms de fichiers sans extensions.') ;
-    end
-end
-
-s1 = [s1 [' ';' ']] ;
-file1 = [s1(1,1:min(find(abs(s1(1,:)) == 32))-1) '.BIN'] ;
-file2 = [s1(2,1:min(find(abs(s1(2,:)) == 32))-1) '.BIN'] ;
-
-fid=fopen(file1,'r') ;
-n1 = fread(fid,1,'int') ;
-n2 = fread(fid,1,'int') ;
-n3 = fread(fid,1,'int') ;
-lag1 = fread(fid,4,'int') ;
-nvx = fread(fid,[n1,n3],'int') ;
-x = fread(fid,[n1,n2],'float64') ;
-fclose(fid) ;
-nvx = char(nvx) ;
-
-fid=fopen(file2,'r') ;
-n1 = fread(fid,1,'int') ;
-n2 = fread(fid,1,'int') ;
-n3 = fread(fid,1,'int') ;
-lag2 = fread(fid,4,'int') ;
-nvy = fread(fid,[n1,n3],'int') ;
-y = fread(fid,[n1,n2],'float64') ;
-fclose(fid) ;
-nvy = char(nvy) ;
-
-if size(x,1) ~= size(y,1)
-    error ('FTEST: The two files don''t have the same number of variables.');
-end
-
-for i = 1:size(x,1)
-    if ~ strcmp(nvx(i,:),nvy(i,:))
-        error ('FTEST: The two files don''t have the same  variables.') ;
-    end
-end
-
-if nnz(lag1 - lag2) > 0
-    error ('FTEST: Leads and lags aren''t the same in both files.') ;
-end
-
-j = zeros(size(s2,1),1);
-for i=1:size(s2,1)
-    k = strmatch(s2(i,:),nvx,'exact') ;
-    if isempty(k)
-        t = ['FTEST: Variable ' s2(i) 'doesn''t exist'] ;
-        error (t) ;
-    else
-        j(i) =k;
-    end
-end
-
-y = y(j,:) ;
-x = x(j,:) ;
-
-%06/18/01 MJ replaced beastr by strmatch
diff --git a/matlab/gcompare.m b/matlab/gcompare.m
deleted file mode 100644
index 3e7644fc916c02cb324ab1888efee3a81cd86405..0000000000000000000000000000000000000000
--- a/matlab/gcompare.m
+++ /dev/null
@@ -1,52 +0,0 @@
-function gcompare(s1,s2)
-% GCOMPARE :    GCOMPARE ( [ 'file1' ; 'file2' ] , [ 'var1' ; 'var2' ...] )
-%               This optional command plots the trajectories of a list of
-%               variables in two different simulations. One plot is drawn
-%               for each variable. The trajectories must have been previously
-%               saved by the instruction DYNASAVE. The simulation in file1
-%               is refered to as the base simulation.
-
-% Copyright © 2001-2017 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 <https://www.gnu.org/licenses/>.
-
-global options_ M_
-global nvx nvy x y lag1
-
-ftest(s1,s2) ;
-
-ix = [1-lag1(1):size(x,2)-lag1(1)]' ;
-i = [lag1(1):size(ix,1)-lag1(2)+1]' ;
-
-if options_.smpl == 0
-    i = [M_.maximum_lag:size(y,2)]' ;
-else
-    i = [options_.smpl(1):options_.smpl(2)]' ;
-end
-
-for k = 1:size(x,1)
-    figure ;
-    plot (ix(i),x(k,i),ix(i),y(k,i)) ;
-    xlabel (['Periods']) ;
-    title (['Variable ' s2(k,:)]) ;
-    l = min(i) + 1;
-    ll = max(i) - 1 ;
-    text (l,x(k,l),s1(1,:)) ;
-    text (ll,y(k,ll),s1(2,:)) ;
-end
-
-% 06/18/01 MJ corrected treatment of options_.smpl
-% 06/24/01 MJ removed color specification
diff --git a/matlab/get_Hessian.m b/matlab/get_Hessian.m
deleted file mode 100644
index e55d79093a100312463dcadffc48103c1d9b55a4..0000000000000000000000000000000000000000
--- a/matlab/get_Hessian.m
+++ /dev/null
@@ -1,255 +0,0 @@
-function [Hess] = get_Hessian(T,R,Q,H,P,Y,DT,DYss,DOm,DH,DP,D2T,D2Yss,D2Om,D2H,D2P,start,mf,kalman_tol,riccati_tol)
-% function [Hess] = get_Hessian(T,R,Q,H,P,Y,DT,DYss,DOm,DH,DP,D2T,D2Yss,D2Om,D2H,D2P,start,mf,kalman_tol,riccati_tol)
-%
-% computes the hessian matrix of the log-likelihood function of
-% a state space model (notation as in kalman_filter.m in DYNARE
-% Thanks to  Nikolai Iskrev
-%
-% NOTE: the derivative matrices (DT,DR ...) are 3-dim. arrays with last
-% dimension equal to the number of structural parameters
-% NOTE: the derivative matrices (D2T,D2Om ...) are 4-dim. arrays with last
-% two dimensions equal to the number of structural parameters
-
-% Copyright © 2011-2017 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 <https://www.gnu.org/licenses/>.
-
-
-k = size(DT,3);                                 % number of structural parameters
-smpl = size(Y,2);                               % Sample size.
-pp   = size(Y,1);                               % Maximum number of observed variables.
-mm   = size(T,2);                               % Number of state variables.
-a    = zeros(mm,1);                             % State vector.
-Om   = R*Q*transpose(R);                        % Variance of R times the vector of structural innovations.
-t    = 0;                                       % Initialization of the time index.
-oldK = 0;
-notsteady   = 1;                                % Steady state flag.
-F_singular  = 1;
-
-Hess  = zeros(k,k);                             % Initialization of the Hessian
-Da    = zeros(mm,k);                             % State vector.
-Dv = zeros(length(mf),k);
-D2a    = zeros(mm,k,k);                             % State vector.
-D2v = zeros(length(mf),k,k);
-
-C = zeros(length(mf),mm);
-for ii=1:length(mf); C(ii,mf(ii))=1;end         % SELECTION MATRIX IN MEASUREMENT EQ. (FOR WHEN IT IS NOT CONSTANT)
-dC = zeros(length(mf),mm,k);
-d2C = zeros(length(mf),mm,k,k);
-
-s   = zeros(pp,1);                      % CONSTANT TERM IN MEASUREMENT EQ. (FOR WHEN IT IS NOT CONSTANT)
-ds  = zeros(pp,1,k);
-d2s = zeros(pp,1,k,k);
-
-%     for ii = 1:k
-%         DOm = DR(:,:,ii)*Q*transpose(R) + R*DQ(:,:,ii)*transpose(R) + R*Q*transpose(DR(:,:,ii));
-%     end
-
-while notsteady & t<smpl
-    t  = t+1;
-    v  = Y(:,t)-a(mf);
-    F  = P(mf,mf) + H;
-    if rcond(F) < kalman_tol
-        if ~all(abs(F(:))<kalman_tol)
-            return
-        else
-            a = T*a;
-            P = T*P*transpose(T)+Om;
-        end
-    else
-        F_singular = 0;
-        iF     = inv(F);
-        K      = P(:,mf)*iF;
-
-        [DK,DF,DP1] = computeDKalman(T,DT,DOm,P,DP,DH,mf,iF,K);
-        [D2K,D2F,D2P1] = computeD2Kalman(T,DT,D2T,D2Om,P,DP,D2P,DH,mf,iF,K,DK);
-        tmp = (a+K*v);
-
-        for ii = 1:k
-            Dv(:,ii)   = -Da(mf,ii) - DYss(mf,ii);
-            %  dai = da(:,:,ii);
-            dKi  = DK(:,:,ii);
-            diFi = -iF*DF(:,:,ii)*iF;
-            dtmpi = Da(:,ii)+dKi*v+K*Dv(:,ii);
-
-
-            for jj = 1:ii
-                dFj    = DF(:,:,jj);
-                diFj   = -iF*DF(:,:,jj)*iF;
-                dKj  = DK(:,:,jj);
-                d2Kij  = D2K(:,:,jj,ii);
-                d2Fij  = D2F(:,:,jj,ii);
-                d2iFij = -diFi*dFj*iF -iF*d2Fij*iF -iF*dFj*diFi;
-                dtmpj = Da(:,jj)+dKj*v+K*Dv(:,jj);
-
-                d2vij  = -D2Yss(mf,jj,ii)  - D2a(mf,jj,ii);
-                d2tmpij = D2a(:,jj,ii) + d2Kij*v + dKj*Dv(:,ii) + dKi*Dv(:,jj) + K*d2vij;
-                D2a(:,jj,ii) = D2T(:,:,jj,ii)*tmp + DT(:,:,jj)*dtmpi + DT(:,:,ii)*dtmpj + T*d2tmpij;
-
-                Hesst(ii,jj) = getHesst_ij(v,Dv(:,ii),Dv(:,jj),d2vij,iF,diFi,diFj,d2iFij,dFj,d2Fij);
-            end
-            Da(:,ii)   = DT(:,:,ii)*tmp + T*dtmpi;
-        end
-        %                     vecDPmf = reshape(DP(mf,mf,:),[],k);
-        %                     iPmf = inv(P(mf,mf));
-        if t>=start
-            Hess = Hess + Hesst;
-        end
-        a      = T*(a+K*v);
-        P      = T*(P-K*P(mf,:))*transpose(T)+Om;
-        DP     = DP1;
-        D2P     = D2P1;
-    end
-    notsteady = max(max(abs(K-oldK))) > riccati_tol;
-    oldK = K;
-end
-
-if F_singular
-    error('The variance of the forecast error remains singular until the end of the sample')
-end
-
-
-if t < smpl
-    t0 = t+1;
-    while t < smpl
-        t = t+1;
-        v = Y(:,t)-a(mf);
-        tmp = (a+K*v);
-        for ii = 1:k
-            Dv(:,ii)   = -Da(mf,ii)-DYss(mf,ii);
-            dKi  = DK(:,:,ii);
-            diFi = -iF*DF(:,:,ii)*iF;
-            dtmpi = Da(:,ii)+dKi*v+K*Dv(:,ii);
-
-            for jj = 1:ii
-                dFj    = DF(:,:,jj);
-                diFj   = -iF*DF(:,:,jj)*iF;
-                dKj  = DK(:,:,jj);
-                d2Kij  = D2K(:,:,jj,ii);
-                d2Fij  = D2F(:,:,jj,ii);
-                d2iFij = -diFi*dFj*iF -iF*d2Fij*iF -iF*dFj*diFi;
-                dtmpj = Da(:,jj)+dKj*v+K*Dv(:,jj);
-
-                d2vij  = -D2Yss(mf,jj,ii)  - D2a(mf,jj,ii);
-                d2tmpij = D2a(:,jj,ii) + d2Kij*v + dKj*Dv(:,ii) + dKi*Dv(:,jj) + K*d2vij;
-                D2a(:,jj,ii) = D2T(:,:,jj,ii)*tmp + DT(:,:,jj)*dtmpi + DT(:,:,ii)*dtmpj + T*d2tmpij;
-
-                Hesst(ii,jj) = getHesst_ij(v,Dv(:,ii),Dv(:,jj),d2vij,iF,diFi,diFj,d2iFij,dFj,d2Fij);
-            end
-            Da(:,ii)   = DT(:,:,ii)*tmp + T*dtmpi;
-        end
-        if t>=start
-            Hess = Hess + Hesst;
-        end
-        a = T*(a+K*v);
-    end
-    %         Hess = Hess + .5*(smpl+t0-1)*(vecDPmf' * kron(iPmf,iPmf) * vecDPmf);
-    %         for ii = 1:k;
-    %             for jj = 1:ii
-    %              H(ii,jj) = trace(iPmf*(.5*DP(mf,mf,ii)*iPmf*DP(mf,mf,jj) + Dv(:,ii)*Dv(:,jj)'));
-    %             end
-    %         end
-end
-
-Hess = Hess + tril(Hess,-1)';
-
-Hess = -Hess/2;
-% end of main function
-
-function Hesst_ij = getHesst_ij(e,dei,dej,d2eij,iS,diSi,diSj,d2iSij,dSj,d2Sij);
-% computes (i,j) term in the Hessian
-
-Hesst_ij = trace(diSi*dSj + iS*d2Sij) + e'*d2iSij*e + 2*(dei'*diSj*e + dei'*iS*dej + e'*diSi*dej + e'*iS*d2eij);
-
-% end of getHesst_ij
-
-function [DK,DF,DP1] = computeDKalman(T,DT,DOm,P,DP,DH,mf,iF,K)
-
-k      = size(DT,3);
-tmp    = P-K*P(mf,:);
-
-for ii = 1:k
-    DF(:,:,ii)  = DP(mf,mf,ii) + DH(:,:,ii);
-    DiF(:,:,ii) = -iF*DF(:,:,ii)*iF;
-    DK(:,:,ii)  = DP(:,mf,ii)*iF + P(:,mf)*DiF(:,:,ii);
-    Dtmp        = DP(:,:,ii) - DK(:,:,ii)*P(mf,:) - K*DP(mf,:,ii);
-    DP1(:,:,ii) = DT(:,:,ii)*tmp*T' + T*Dtmp*T' + T*tmp*DT(:,:,ii)' + DOm(:,:,ii);
-end
-
-% end of computeDKalman
-
-function [d2K,d2S,d2P1] = computeD2Kalman(A,dA,d2A,d2Om,P0,dP0,d2P0,DH,mf,iF,K0,dK0)
-% computes the second derivatives of the Kalman matrices
-% note: A=T in main func.
-
-k      = size(dA,3);
-tmp    = P0-K0*P0(mf,:);
-[ns,no] = size(K0);
-
-% CPC = C*P0*C'; CPC = .5*(CPC+CPC');iF = inv(CPC);
-% APC = A*P0*C';
-% APA = A*P0*A';
-
-
-d2K  = zeros(ns,no,k,k);
-d2S  = zeros(no,no,k,k);
-d2P1 = zeros(ns,ns,k,k);
-
-for ii = 1:k
-    dAi = dA(:,:,ii);
-    dFi = dP0(mf,mf,ii);
-    d2Omi = d2Om(:,:,ii);
-    diFi = -iF*dFi*iF;
-    dKi = dK0(:,:,ii);
-    for jj = 1:k
-        dAj = dA(:,:,jj);
-        dFj = dP0(mf,mf,jj);
-        d2Omj = d2Om(:,:,jj);
-        dFj = dP0(mf,mf,jj);
-        diFj = -iF*dFj*iF;
-        dKj = dK0(:,:,jj);
-
-        d2Aij = d2A(:,:,jj,ii);
-        d2Pij = d2P0(:,:,jj,ii);
-        d2Omij = d2Om(:,:,jj,ii);
-
-        % second order
-
-        d2Fij = d2Pij(mf,mf) ;
-
-        %     d2APC = d2Aij*P0*C' + A*d2Pij*C' + A*P0*d2Cij' + dAi*dPj*C' + dAj*dPi*C' + A*dPj*dCi' + A*dPi*dCj' + dAi*P0*dCj' + dAj*P0*dCi';
-        d2APC = d2Pij(:,mf);
-
-        d2iF = -diFi*dFj*iF -iF*d2Fij*iF -iF*dFj*diFi;
-
-        d2Kij= d2Pij(:,mf)*iF + P0(:,mf)*d2iF + dP0(:,mf,jj)*diFi + dP0(:,mf,ii)*diFj;
-
-        d2KCP = d2Kij*P0(mf,:) + K0*d2Pij(mf,:) + dKi*dP0(mf,:,jj) + dKj*dP0(mf,:,ii) ;
-
-        dtmpi        = dP0(:,:,ii) - dK0(:,:,ii)*P0(mf,:) - K0*dP0(mf,:,ii);
-        dtmpj        = dP0(:,:,jj) - dK0(:,:,jj)*P0(mf,:) - K0*dP0(mf,:,jj);
-        d2tmp = d2Pij - d2KCP;
-
-        d2AtmpA = d2Aij*tmp*A' + A*d2tmp*A' + A*tmp*d2Aij' + dAi*dtmpj*A' + dAj*dtmpi*A' + A*dtmpj*dAi' + A*dtmpi*dAj' + dAi*tmp*dAj' + dAj*tmp*dAi';
-
-        d2K(:,:,ii,jj)  = d2Kij; %#ok<NASGU>
-        d2P1(:,:,ii,jj) = d2AtmpA  + d2Omij;  %#ok<*NASGU>
-        d2S(:,:,ii,jj)  = d2Fij;
-        %     d2iS(:,:,ii,jj) = d2iF;
-    end
-end
-
-% end of computeD2Kalman
diff --git a/matlab/mcpath_function.m b/matlab/mcpath_function.m
deleted file mode 100644
index dc985b382c4dc73c8d09ae0c686211ed1394effe..0000000000000000000000000000000000000000
--- a/matlab/mcpath_function.m
+++ /dev/null
@@ -1,8 +0,0 @@
-function [res,jac,domerr] = mcpath_function(func,z,jacflag,varargin)
-domerr = 0;
-
-if jacflag
-    [res,jac] = func(z,varargin{:});
-else
-    res = func(z,varargin{:});
-end
diff --git a/matlab/metropolis_run_analysis.m b/matlab/metropolis_run_analysis.m
deleted file mode 100644
index 0baa03ab202eee9410655ad7c0a567a420674632..0000000000000000000000000000000000000000
--- a/matlab/metropolis_run_analysis.m
+++ /dev/null
@@ -1,55 +0,0 @@
-function metropolis_run_analysis(M_,basetopt,j)
-%function metropolis_run_analysis(M_,basetopt,j
-% analizes Metropolis runs
-%
-% INPUTS
-%   M_:        (struct)  Model structure
-%   basetopt:  (struct)  Estimated parameters structure
-%   j:         (int)     Index of estimated paramter
-%
-% OUTPUTS
-%   none
-%
-% SPECIAL REQUIREMENTS
-%   none
-
-% Copyright © 2003-2017 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 <https://www.gnu.org/licenses/>.
-
-load([M_.fname '/metropolis/' M_.fname '_mh_history'])
-nblck = record.Nblck;
-ndraws = sum(record.MhDraws(:,1));
-
-logPost = [];
-params = [];
-blck = 1;
-for i=1:record.LastFileNumber
-    fname = [M_.fname '/metropolis/' M_.fname '_mh' int2str(i) '_blck' ...
-             int2str(blck) '.mat'];
-    if exist(fname,'file')
-        o=load(fname);
-        logPost = [logPost; o.logpo2];
-        params  = [params; o.x2];
-    end
-end
-
-figure;
-subplot(2,1,1)
-plot(logPost)
-subplot(2,1,2)
-plot(params(:,j))
-title(['parameter ' int2str(j)])
\ No newline at end of file
diff --git a/matlab/score.m b/matlab/score.m
deleted file mode 100644
index 0ae59bf652f73d68a867b77caf3020e8f6b82bdb..0000000000000000000000000000000000000000
--- a/matlab/score.m
+++ /dev/null
@@ -1,123 +0,0 @@
-function [DLIK] = score(T,R,Q,H,P,Y,DT,DYss,DOm,DH,DP,start,mf,kalman_tol,riccati_tol)
-% function [DLIK] = score(T,R,Q,H,P,Y,DT,DYss,DOm,DH,DP,start,mf,kalman_tol,riccati_tol)
-%
-% computes the derivative of the log-likelihood function of
-% a state space model (notation as in kalman_filter.m in DYNARE
-% thanks to Nikolai Iskrev
-%
-% NOTE: the derivative matrices (DT,DR ...) are 3-dim. arrays with last
-% dimension equal to the number of structural parameters
-
-% Copyright © 2009-2017 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/licen
-
-k = size(DT,3);                                 % number of structural parameters
-smpl = size(Y,2);                               % Sample size.
-mm   = size(T,2);                               % Number of state variables.
-a    = zeros(mm,1);                             % State vector.
-Om   = R*Q*transpose(R);                        % Variance of R times the vector of structural innovations.
-t    = 0;                                       % Initialization of the time index.
-oldK = 0;
-notsteady   = 1;                                % Steady state flag.
-F_singular  = 1;
-
-DLIK  = zeros(k,1);                             % Initialization of the score.
-Da    = zeros(mm,k);                            % State vector.
-Dv    = zeros(length(mf),k);                    % observation vector.
-
-%     for ii = 1:k
-%         DOm = DR(:,:,ii)*Q*transpose(R) + R*DQ(:,:,ii)*transpose(R) + R*Q*transpose(DR(:,:,ii));
-%     end
-
-while notsteady & t<smpl
-    t  = t+1;
-    v  = Y(:,t)-a(mf);
-    F  = P(mf,mf) + H;
-    if rcond(F) < kalman_tol
-        if ~all(abs(F(:))<kalman_tol)
-            return
-        else
-            a = T*a;
-            P = T*P*transpose(T)+Om;
-        end
-    else
-        F_singular = 0;
-        iF     = inv(F);
-        K      = P(:,mf)*iF;
-
-        [DK,DF,DP1] = computeDKalman(T,DT,DOm,P,DP,DH,mf,iF,K);
-        for ii = 1:k
-            Dv(:,ii)   = -Da(mf,ii)-DYss(mf,ii);
-            Da(:,ii)   = DT(:,:,ii)*(a+K*v) + T*(Da(:,ii)+DK(:,:,ii)*v + K*Dv(:,ii));
-            if t>=start
-                DLIK(ii,1)  = DLIK(ii,1) + trace( iF*DF(:,:,ii) ) + 2*Dv(:,ii)'*iF*v - v'*(iF*DF(:,:,ii)*iF)*v;
-            end
-        end
-        a      = T*(a+K*v);
-        P      = T*(P-K*P(mf,:))*transpose(T)+Om;
-        DP     = DP1;
-    end
-    notsteady = max(max(abs(K-oldK))) > riccati_tol;
-    oldK = K;
-end
-
-if F_singular
-    error('The variance of the forecast error remains singular until the end of the sample')
-end
-
-for ii = 1:k
-    tmp0(:,:,ii) = iF*DF(:,:,ii)*iF;
-end
-
-if t < smpl
-    t0 = t+1;
-    while t < smpl
-        t = t+1;
-        v = Y(:,t)-a(mf);
-        for ii = 1:k
-            Dv(:,ii)   = -Da(mf,ii)-DYss(mf,ii);
-            Da(:,ii)   = DT(:,:,ii)*(a+K*v) + T*(Da(:,ii)+DK(:,:,ii)*v + K*Dv(:,ii));
-            if t>=start
-                DLIK(ii,1)  = DLIK(ii,1) + trace( iF*DF(:,:,ii) ) + 2*Dv(:,ii)'*iF*v - v'*(iF*DF(:,:,ii)*iF)*v;
-            end
-        end
-        a = T*(a+K*v);
-    end
-    for ii = 1:k
-        %             DLIK(ii,1)  = DLIK(ii,1) + (smpl-t0+1)*trace( iF*DF(:,:,ii) );
-    end
-
-end
-
-DLIK = DLIK/2;
-
-% end of main function
-
-function [DK,DF,DP1] = computeDKalman(T,DT,DOm,P,DP,DH,mf,iF,K)
-
-k      = size(DT,3);
-tmp    = P-K*P(mf,:);
-
-for ii = 1:k
-    DF(:,:,ii)  = DP(mf,mf,ii) + DH(:,:,ii);
-    DiF(:,:,ii) = -iF*DF(:,:,ii)*iF;
-    DK(:,:,ii)  = DP(:,mf,ii)*iF + P(:,mf)*DiF(:,:,ii);
-    Dtmp        = DP(:,:,ii) - DK(:,:,ii)*P(mf,:) - K*DP(mf,:,ii);
-    DP1(:,:,ii) = DT(:,:,ii)*tmp*T' + T*Dtmp*T' + T*tmp*DT(:,:,ii)' + DOm(:,:,ii);
-end
-
-% end of computeDKalman
diff --git a/matlab/shiftS.m b/matlab/shiftS.m
deleted file mode 100644
index f0f7ba8abff1a37aeee09afc63a5f089793be16a..0000000000000000000000000000000000000000
--- a/matlab/shiftS.m
+++ /dev/null
@@ -1,28 +0,0 @@
-function S = shiftS(S,n)
-%function S = shiftS(S,n)
-%
-% Removes the first n elements of a one dimensional cell array.
-
-% Copyright © 2013-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 <https://www.gnu.org/licenses/>.
-
-if length(S) >= n+1
-    S = S(n+1:end);
-else
-    S = {};
-end
-end
\ No newline at end of file