diff --git a/matlab/dynare_config.m b/matlab/dynare_config.m
index 0769f9e3e..0ed36dde4 100644
--- a/matlab/dynare_config.m
+++ b/matlab/dynare_config.m
@@ -38,6 +38,8 @@ end
 dynareroot = strrep(which('dynare.m'),'dynare.m','');
 
 addpath([dynareroot '/distributions/'])
+addpath([dynareroot '/kalman/'])
+
 % Add path to distribution-related function if under Matlab
 % without the statistics toolbox
 if ~exist('OCTAVE_VERSION') && isempty(ver('stats'))
diff --git a/matlab/global_initialization.m b/matlab/global_initialization.m
index d3489b0b0..38adfd6d2 100644
--- a/matlab/global_initialization.m
+++ b/matlab/global_initialization.m
@@ -145,6 +145,7 @@ function global_initialization()
   options_.first_obs = 1;
   options_.kalman_algo = 0;
   options_.kalman_tol = 1e-12;
+  options_.riccati_tol = 1e-6;
   options_.lik_algo = 1;
   options_.lik_init = 1;
   options_.load_mh_file = 0;
diff --git a/matlab/kalman/build_selection_matrix.m b/matlab/kalman/build_selection_matrix.m
new file mode 100644
index 000000000..bd47b72c1
--- /dev/null
+++ b/matlab/kalman/build_selection_matrix.m
@@ -0,0 +1,24 @@
+function Z = build_selection_matrix(mf,m,p)
+% Builds the selection matrix of the measurement equation from the vector
+% of indices mf.
+
+% Copyright (C) 2004-2008 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/>.    
+    Z = zeros(p,m);
+    for i=1:p
+        Z(i,mf(i)) = 1;
+    end
\ No newline at end of file
diff --git a/matlab/kalman/kalman_filter.m b/matlab/kalman/kalman_filter.m
new file mode 100644
index 000000000..c4e926446
--- /dev/null
+++ b/matlab/kalman/kalman_filter.m
@@ -0,0 +1,96 @@
+function [LIK, lik] = kalman_filter(T,R,Q,H,P,Y,start,mf,kalman_tol,riccati_tol)
+% Computes the likelihood of a stationnary state space model without measurement error, in the case of a non-singular var-cov matrix
+%
+% INPUTS
+%    T                      [double]    mm*mm transition matrix of the state equation.
+%    R                      [double]    mm*rr matrix, mapping structural innovations to state variables.
+%    Q                      [double]    rr*rr covariance matrix of the structural innovations.
+%    H                      [double]    pp*pp (or 1*1 =0 if no measurement error) covariance matrix of the measurement errors. 
+%    P                      [double]    mm*mm variance-covariance matrix with stationary variables
+%    Y                      [double]    pp*smpl matrix of (detrended) data, where pp is the maximum number of observed variables.
+%    start                  [integer]   scalar, likelihood evaluation starts at 'start'.
+%    mf                     [integer]   pp*1 vector of indices.
+%    kalman_tol             [double]    scalar, tolerance parameter (rcond).
+%    riccati_tol            [double]    scalar, tolerance parameter (riccati iteration).
+%
+% OUTPUTS
+%    LIK        [double]    scalar, likelihood
+%    lik        [double]    vector, density of observations in each period.
+%
+% REFERENCES
+%   See "Filtering and Smoothing of State Vector for Diffuse State Space
+%   Models", S.J. Koopman and J. Durbin (2003, in Journal of Time Series
+%   Analysis, vol. 24(1), pp. 85-98).
+%
+% NOTES
+%   The vector "lik" is used to evaluate the jacobian of the likelihood.
+
+% Copyright (C) 2004-2008 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/>.
+
+    smpl = size(Y,2);                               % Sample size.
+    mm   = size(T,2);                               % Number of state variables.
+    pp   = size(Y,1);                               % Maximum number of observed variables.
+    a    = zeros(mm,1);                             % State vector.
+    dF   = 1;                                       % det(F).
+    QQ   = R*Q*transpose(R);                        % Variance of R times the vector of structural innovations.
+    t    = 0;                                       % Initialization of the time index.
+    lik  = zeros(smpl+1,1);                         % Initialization of the vector gathering the densities.
+    LIK  = Inf;                                     % Default value of the log likelihood.
+    oldK = 0;
+    number_of_observations = smpl*pp;
+    lik(smpl+1) = number_of_observations*log(2*pi); % Log likelihood constant.
+    notsteady   = 1;                                % Steady state flag.
+    F_singular  = 1;
+
+    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)+QQ;
+            end
+        else
+            F_singular = 0;
+            dF     = det(F);
+            iF     = inv(F);
+            lik(t) = log(dF)+transpose(v)*iF*v;
+            K      = P(:,mf)*iF;
+            a      = T*(a+K*v);
+            P      = T*(P-K*P(mf,:))*transpose(T)+QQ;
+        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
+
+    reste = smpl-t;
+    while t < smpl
+        t = t+1;
+        v = Y(:,t)-a(mf);
+        a = T*(a+K*v);
+        lik(t) = transpose(v)*iF*v;
+    end
+    lik(t) = lik(t) + reste*log(dF);
+    LIK    = .5*(sum(lik(start:end))-(start-1)*lik(smpl+1)/smpl);% Minus the log-likelihood.
\ No newline at end of file
diff --git a/matlab/kalman/missing_observations_kalman_filter.m b/matlab/kalman/missing_observations_kalman_filter.m
new file mode 100644
index 000000000..17bd02efa
--- /dev/null
+++ b/matlab/kalman/missing_observations_kalman_filter.m
@@ -0,0 +1,110 @@
+function  [LIK, lik] = missing_observations_kalman_filter(T,R,Q,H,P,Y,start,mf,kalman_tol,riccati_tol,data_index,number_of_observations,no_more_missing_observations)
+% Computes the likelihood of a state space model in the case with missing observations.
+%
+% INPUTS
+%    T                      [double]    mm*mm transition matrix of the state equation.
+%    R                      [double]    mm*rr matrix, mapping structural innovations to state variables.
+%    Q                      [double]    rr*rr covariance matrix of the structural innovations.
+%    H                      [double]    pp*pp (or 1*1 =0 if no measurement error) covariance matrix of the measurement errors.    
+%    P                      [double]    mm*mm variance-covariance matrix of the initial state vector.
+%    Y                      [double]    pp*smpl matrix of detrended data, where pp is the maximum number of observed variables
+%    trend                  [double]    pp*smpl matrix
+%    start                  [integer]   scalar, likelihood evaluation starts at 'start'.
+%    mf                     [integer]   pp*1 vector of indices.
+%    kalman_tol             [double]    scalar, tolerance parameter (rcond).
+%    riccati_tol            [double]    scalar, tolerance parameter (riccati iteration).
+%    data_index             [cell]      1*smpl cell of column vectors of indices.
+%    number_of_observations [integer]   saclar.
+%    
+% OUTPUTS
+%    LIK        [double]    scalar, likelihood
+%    lik        [double]    vector, density of observations in each period.
+%        
+% REFERENCES
+%   See "Filtering and Smoothing of State Vector for Diffuse State Space
+%   Models", S.J. Koopman and J. Durbin (2003, in Journal of Time Series 
+%   Analysis, vol. 24(1), pp. 85-98). 
+%
+% NOTES
+%   The vector "lik" is used to evaluate the jacobian of the likelihood.   
+    
+% Copyright (C) 2004-2008 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/>.
+    
+    smpl = size(Y,2);                               % Sample size.
+    mm   = size(T,2);                               % Number of state variables.
+    pp   = size(Y,1);                               % Maximum number of observed variables.
+    a    = zeros(mm,1);                             % State vector.       
+    dF   = 1;                                       % det(F).
+    QQ   = R*Q*transpose(R);                        % Variance of R times the vector of structural innovations. 
+    t    = 0;                                       % Initialization of the time index.
+    lik  = zeros(smpl+1,1);                         % Initialization of the vector gathering the densities.
+    LIK  = Inf;                                     % Default value of the log likelihood.
+    oldK = 0;
+    lik(smpl+1) = number_of_observations*log(2*pi); % Log likelihood constant.  
+    notsteady   = 1;                                % Steady state flag.
+    F_singular  = 1;
+    
+    while notsteady & t<smpl 
+        t  = t+1;
+        if isempty(data_index{t})
+            a = T*a;
+            P = T*P*transpose(T)+QQ;
+        else
+            MF = mf(data_index{t});                 % Set the selection for observed variables. 
+            v  = Y(data_index{t},t)-a(MF);
+            if ~isscalar(H)                         % => Errors in the measurement equation.
+                F = P(MF,MF) + H(data_index{t},data_index{t}); 
+            else% => 
+                % case 1. No errors in the measurement (H=0) and more than one variable is observed in this state space model. 
+                % case 2. Errors in the measurement equation, but only one variable is observed in this state-space model.
+                F = P(MF,MF)+H;
+            end
+            if rcond(F) < kalman_tol
+                if ~all(abs(F(:))<kalman_tol)
+                    return
+                else
+                    a = T*a;
+                    P = T*P*transpose(T)+QQ;
+                end
+            else
+                F_singular = 0;
+                dF     = det(F);
+                iF     = inv(F);
+                lik(t) = log(dF)+transpose(v)*iF*v;
+                K      = P(:,MF)*iF;
+                a      = T*(a+K*v);
+                P      = T*(P-K*P(MF,:))*transpose(T)+QQ;
+            end
+            notsteady = ~( (t>no_more_missing_observations) && (max(max(abs(K-oldK)))<riccati_tol) );
+            oldK = K;
+        end
+    end
+    
+    if F_singular
+        error('The variance of the forecast error remains singular until the end of the sample')
+    end
+    
+    reste = smpl-t;
+    while t < smpl
+        t = t+1;
+        v = Y(:,t)-a(mf);
+        a = T*(a+K*v);
+        lik(t) = transpose(v)*iF*v;
+    end
+    lik(t) = lik(t) + reste*log(dF);
+    LIK    = .5*(sum(lik(start:end))-(start-1)*lik(smpl+1)/smpl);% Minus the log-likelihood.
\ No newline at end of file
diff --git a/matlab/kalman/steady_state_kalman_gain.m b/matlab/kalman/steady_state_kalman_gain.m
new file mode 100644
index 000000000..2bede6e0c
--- /dev/null
+++ b/matlab/kalman/steady_state_kalman_gain.m
@@ -0,0 +1,50 @@
+function [K,iF,P] = steady_state_kalman_gain(T,R,Q,H,mf)
+% Given the invariant state space representation of a model, this
+% function computes the gain matrix and the covariance matrix of the
+% state vector at the steady state of the kalman filter. 
+% 
+% INPUTS 
+%   T   [double]    m*m transition matrix of the state vector.  
+%   R   [double]    m*q matrix (q is the number of structural innovations).
+%   Q   [double]    q*q covariance matrix of the structural innovations.
+%   H   [double]    p*p covariance matrix of the measurement error.
+%   mf  [integer]   p*1 vector, indices for the observed variables
+%    
+% OUTPUTS 
+%   K   [double]    kalman gain matrix.
+%   P   [double]    covariance matrix of the state vector.
+%               
+% SPECIAL REQUIREMENTS
+%   Needs a solver for Riccati equations (dare.m)
+
+% Copyright (C) 2004-2008 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/>.
+    
+m = length(T);
+p = length(mf);
+Z = build_selection_matrix(mf,m,p);
+
+if isempty(H)
+    H = zeros(p,p);
+end
+
+QQ = R*Q*transpose(R);
+
+P  = dare(T,transpose(Z),QQ,H);
+
+iF = inv(Z*P*transpose(Z)+H);
+K  = T*P*transpose(Z)*iF;
\ No newline at end of file