diff --git a/matlab/DsgeSmoother.m b/matlab/DsgeSmoother.m
index 4f37af83231e4cb07808ee7c6fd8d80442998d94..8cbe930bd6de631c2d541d2f9c54bccbda9c51c5 100644
--- a/matlab/DsgeSmoother.m
+++ b/matlab/DsgeSmoother.m
@@ -1,4 +1,4 @@
-function [alphahat,etahat,epsilonhat,ahat,SteadyState,trend_coeff] = DsgeSmoother(xparam1,gend,Y)
+function [alphahat,etahat,epsilonhat,ahat,SteadyState,trend_coeff,aK] = DsgeSmoother(xparam1,gend,Y)
 % stephane.adjemian@cepremap.cnrs.fr [09-07-2004]
 %
 % Adapted from mj_optmumlik.m
@@ -10,41 +10,8 @@ function [alphahat,etahat,epsilonhat,ahat,SteadyState,trend_coeff] = DsgeSmoothe
   nobs 		= size(options_.varobs,1);
   smpl        = size(Y,2);
 
-  Q = M_.Sigma_e;
-  for i=1:estim_params_.nvx
-    k =estim_params_.var_exo(i,1);
-    Q(k,k) = xparam1(i)*xparam1(i);
-  end
-  offset = estim_params_.nvx;
-  if estim_params_.nvn
-    H = zeros(nobs,nobs);
-    for i=1:estim_params_.nvn
-      k = estim_params_.var_endo(i,1);
-      H(k,k) = xparam1(i+offset)*xparam1(i+offset);
-    end
-  end	
-  offset = offset+estim_params_.nvn;
-  for i=1:estim_params_.ncx
-    k1 =estim_params_.corrx(i,1);
-    k2 =estim_params_.corrx(i,2);
-    Q(k1,k2) = xparam1(i+offset)*sqrt(Q(k1,k1)*Q(k2,k2));
-    Q(k2,k1) = Q(k1,k2);
-  end
-  offset = offset+estim_params_.ncx;
+  set_all_parameters(xparam1);
 
-  if estim_params_.nvn & estim_params_.ncn 
-    for i=1:estim_params_.ncn
-      k1 = options_.lgyidx2varobs(estim_params_.corrn(i,1));
-      k2 = options_.lgyidx2varobs(estim_params_.corrn(i,2));
-      H(k1,k2) = xparam1(i+offset)*sqrt(H(k1,k1)*H(k2,k2));
-      H(k2,k1) = H(k1,k2);
-    end
-    offset = offset+estim_params_.ncn;
-  end	
-  for i=1:estim_params_.np
-    M_.params(estim_params_.param_vals(i,1)) = xparam1(i+offset);
-  end
-  M_.Sigma_e = Q;
   %------------------------------------------------------------------------------
   % 2. call model setup & reduction program
   %------------------------------------------------------------------------------
@@ -76,6 +43,9 @@ function [alphahat,etahat,epsilonhat,ahat,SteadyState,trend_coeff] = DsgeSmoothe
   %  alors il suffit de poser Pstar comme la solution de l'�uation de Lyapounov et
   %  Pinf=[].
   %
+  Q = M_.Sigma_e;
+  H = M_.H;
+  
   if options_.lik_init == 1		% Kalman filter
     Pstar = lyapunov_symm(T,R*Q*transpose(R));
     Pinf	= [];
@@ -101,20 +71,20 @@ function [alphahat,etahat,epsilonhat,ahat,SteadyState,trend_coeff] = DsgeSmoothe
   % -----------------------------------------------------------------------------
   if estim_params_.nvn
     if options_.kalman_algo == 1
-      [alphahat,epsilonhat,etahat,ahat] = DiffuseKalmanSmootherH1(T,R,Q,H,Pinf,Pstar,Y,trend,nobs,np,smpl,mf);
+      [alphahat,epsilonhat,etahat,ahat,aK] = DiffuseKalmanSmootherH1(T,R,Q,H,Pinf,Pstar,Y,trend,nobs,np,smpl,mf);
       if all(alphahat(:)==0)
-	[alphahat,epsilonhat,etahat,ahat] = DiffuseKalmanSmootherH3(T,R,Q,H,Pinf,Pstar,Y,trend,nobs,np,smpl,mf);
+	[alphahat,epsilonhat,etahat,ahat,aK] = DiffuseKalmanSmootherH3(T,R,Q,H,Pinf,Pstar,Y,trend,nobs,np,smpl,mf);
       end
     elseif options_.kalman_algo == 3
-      [alphahat,epsilonhat,etahat,ahat] = DiffuseKalmanSmootherH3(T,R,Q,H,Pinf,Pstar,Y,trend,nobs,np,smpl,mf);
+      [alphahat,epsilonhat,etahat,ahat,aK] = DiffuseKalmanSmootherH3(T,R,Q,H,Pinf,Pstar,Y,trend,nobs,np,smpl,mf);
     end
   else
     if options_.kalman_algo == 1
-      [alphahat,etahat,ahat] = DiffuseKalmanSmoother1(T,R,Q,Pinf,Pstar,Y,trend,nobs,np,smpl,mf);
+      [alphahat,etahat,ahat,aK] = DiffuseKalmanSmoother1(T,R,Q,Pinf,Pstar,Y,trend,nobs,np,smpl,mf);
       if all(alphahat(:)==0)
 	[alphahat,etahat,ahat] = DiffuseKalmanSmoother3(T,R,Q,Pinf,Pstar,Y,trend,nobs,np,smpl,mf);
       end
     elseif options_.kalman_algo == 3
-      [alphahat,etahat,ahat] = DiffuseKalmanSmoother3(T,R,Q,Pinf,Pstar,Y,trend,nobs,np,smpl,mf);
+      [alphahat,etahat,ahat,aK] = DiffuseKalmanSmoother3(T,R,Q,Pinf,Pstar,Y,trend,nobs,np,smpl,mf);
     end
   end
diff --git a/matlab/GetOneDraw.m b/matlab/GetOneDraw.m
index 20df9a9f2fb53b237657e9e06f48562422c700e1..92d72e75aa04485c2dd5b9ce2b1316a85b7b3ccf 100644
--- a/matlab/GetOneDraw.m
+++ b/matlab/GetOneDraw.m
@@ -1,18 +1,9 @@
-function xparams = GetOneDraw(NumberOfDraws,FirstMhFile,LastMhFile,FirstLine,nn,DirectoryName)
+function xparams = GetOneDraw(type)
 % stephane.adjemian@ens.fr [09-25-2005]
-global options_ M_
 
-ChainNumber = ceil(rand*options_.mh_nblck);
-DrawNumber  = ceil(rand*NumberOfDraws);
-
-if DrawNumber <= nn-FirstLine+1
-  MhFilNumber = FirstMhFile;
-  MhLine = FirstLine+DrawNumber-1;
-else
-  DrawNumber  = DrawNumber-(nn-FirstLine+1);
-  MhFilNumber = FirstMhFile+ceil(DrawNumber/nn); 
-  MhLine = DrawNumber-(MhFilNumber-FirstMhFile-1)*nn;
-end
-
-load([DirectoryName '/' M_.fname '_mh' int2str(MhFilNumber) '_blck' int2str(ChainNumber) '.mat' ],'x2');
-xparams = x2(MhLine,:);
\ No newline at end of file
+  switch type
+   case 'posterior'
+    xparams = metropolis_draw(0);
+   case 'prior'
+    xparams = prior_draw(0);
+  end
\ No newline at end of file
diff --git a/matlab/PosteriorIRF.m b/matlab/PosteriorIRF.m
index 557dda81fd7129bc95fc77cd1821e3277a544c22..43a57e8bacf8efe415bd290d3871d7b4aa41b852 100644
--- a/matlab/PosteriorIRF.m
+++ b/matlab/PosteriorIRF.m
@@ -14,16 +14,6 @@ MaxNumberOfPlotPerFigure = 4;% The square root must be an integer!
 nn = sqrt(MaxNumberOfPlotPerFigure);
 %%
 CheckPath('Plots/IRFs');
-CheckPath('metropolis/IRFs');
-DirectoryName = CheckPath('metropolis');
-load([ DirectoryName '/'  M_.fname '_mh_history'])
-FirstMhFile = record.KeepedDraws.FirstMhFile;
-FirstLine = record.KeepedDraws.FirstLine; 
-TotalNumberOfMhFiles = sum(record.MhDraws(:,2)); LastMhFile = TotalNumberOfMhFiles; 
-TotalNumberOfMhDraws = sum(record.MhDraws(:,1));
-NumberOfDraws = TotalNumberOfMhDraws-floor(options_.mh_drop*TotalNumberOfMhDraws);
-clear record;
-MAX_nruns = ceil(options_.MaxNumberOfBytes/(npar+2)/8);
 MAX_nirfs = ceil(options_.MaxNumberOfBytes/(options_.irf*length(oo_.steady_state)*M_.exo_nbr)/8)+50;
 %%
 B = round(0.25*NumberOfDraws);
@@ -39,26 +29,9 @@ end
 for b=1:B
   irun = irun+1;
   deep = GetOneDraw(NumberOfDraws,FirstMhFile,LastMhFile,FirstLine,MAX_nruns,DirectoryName);
-  M_.params(estim_params_.param_vals(:,1)) = deep(offset+1:end);
+  set_parameters(deep);
   dr = resol(oo_.steady_state,0);
-  if nvx
-    ip = 1;
-    for i=1:nvx
-      k = estim_params_.var_exo(i,1);
-      M_.Sigma_e(k,k) = deep(ip)*deep(ip);
-      ip = ip+1;
-    end
-  end
-  if ncx
-    ip = nvx+nvn+1;
-    for i=1:ncx
-      k1 = estim_params_.corrx(i,1);
-      k2 = estim_params_.corrx(i,2);
-      M_.Sigma_e(k1,k2) = deep(ip)*sqrt(M_.Sigma_e(k1,k1)*M_.Sigma_e(k2,k2));
-      M_.Sigma_e(k2,k1) = M_.Sigma_e(k1,k2);
-      ip = ip+1;
-    end
-  end
+
   SS(M_.exo_names_orig_ord,M_.exo_names_orig_ord) = M_.Sigma_e+1e-14*eye(M_.exo_nbr);
   SS = transpose(chol(SS));
   for i = 1:M_.exo_nbr
diff --git a/matlab/PosteriorSmoother.m b/matlab/PosteriorSmoother.m
new file mode 100644
index 0000000000000000000000000000000000000000..34812486757bf5bd07ca10d389e92be7462f975c
--- /dev/null
+++ b/matlab/PosteriorSmoother.m
@@ -0,0 +1,147 @@
+function PosteriorSmoother(Y,gend)
+% stephane.adjemian@ens.fr [09-25-2005]
+global options_ estim_params_ oo_ M_
+
+nvx  = estim_params_.nvx;
+nvn  = estim_params_.nvn;
+ncx  = estim_params_.ncx;
+ncn  = estim_params_.ncn;
+np   = estim_params_.np ;
+npar = nvx+nvn+ncx+ncn+np;
+offset = npar-np;
+naK = length(options_.filter_step_ahead);
+%%
+MaxNumberOfPlotPerFigure = 4;% The square root must be an integer!
+nn = sqrt(MaxNumberOfPlotPerFigure);
+%%
+CheckPath('Plots/');
+DirectoryName = CheckPath('metropolis');
+load([ DirectoryName '/'  M_.fname '_mh_history'])
+FirstMhFile = record.KeepedDraws.FirstMhFile;
+FirstLine = record.KeepedDraws.FirstLine; 
+TotalNumberOfMhFiles = sum(record.MhDraws(:,2)); LastMhFile = TotalNumberOfMhFiles; 
+TotalNumberOfMhDraws = sum(record.MhDraws(:,1));
+NumberOfDraws = TotalNumberOfMhDraws-floor(options_.mh_drop*TotalNumberOfMhDraws);
+clear record;
+MAX_nruns = ceil(options_.MaxNumberOfBytes/(npar+2)/8);
+MAX_nsmoo = ceil(MaxNumberOfBytes/((endo_nbr)*gend)/8);
+MAX_ninno = ceil(MaxNumberOfBytes/(exo_nbr*gend)/8);
+MAX_nerro = ceil(MaxNumberOfBytes/(size(options_.varobs,1)*gend)/8);
+MAX_naK   = ceil(MaxNumberOfBytes/(size(options_.varobs,1)*length(options_.filter_step_ahead)*gend)/8);
+%%
+B = round(0.25*NumberOfDraws);
+%%
+varlist = options_.varlist;
+if isempty(varlist)
+  varlist = M_.endo_names;
+  SelecVariables = transpose(1:M_.endo_nbr);
+  nvar = M_.endo_nbr;
+else
+  nvar = size(varlist,1);
+  SelecVariables = [];
+  for i=1:nvar
+    if ~isempty(strmatch(varlist(i,:),M_.endo_names,'exact'))
+      SelecVariables = [SelecVariables;strmatch(varlist(i,:),M_.endo_names,'exact')];
+    end
+  end
+end
+
+irun1 = 1;
+irun2 = 1;
+irun3 = 1;
+irun4 = 1;
+ifil1 = 1;
+ifil2 = 1;
+ifil3 = 1;
+ifil4 = 1;
+h = waitbar(0,'Bayesian smoother...');
+if B >= MAX_nirfs 
+  stock_irf = zeros(options_.irf,M_.endo_nbr,M_.exo_nbr,MAX_nirfs);
+else
+  stock_irf = zeros(options_.irf,M_.endo_nbr,M_.exo_nbr,B);
+end
+if options_.smoother
+  if B <= MAX_nsmoo
+    stock_smooth = zeros(endo_nbr,gend,B);
+  else
+    stock_smooth = zeros(endo_nbr,gend,MAX_nsmoo);
+  end
+  if B <= MAX_ninno 
+    stock_innov  = zeros(exo_nbr,gend,B);
+  else
+    stock_innov  = zeros(exo_nbr,gend,MAX_ninno);
+  end
+  if nvn & B <= MAX_nerro
+    stock_error = zeros(nvobs,gend,B);
+  else nvn & B > MAX_nerro
+    stock_error = zeros(nvobs,gend,MAX_nerro);
+  end
+end
+if options_.filter_step_ahead ~= 0
+  if B <= MAX_naK
+    stock_filter = zeros(naK,endo_nbr,gend+1,B);
+  else
+    stock_filter = zeros(naK,endo_nbr,gend+1,MAX_naK);
+  end
+end
+for b=1:B
+  deep = GetOneDraw(NumberOfDraws,FirstMhFile,LastMhFile,FirstLine,MAX_nruns,DirectoryName);
+  set_all_parameters(deep);
+  dr = resol(oo_.steady_state,0);
+  [alphahat,etahat,epsilonhat,ahat,SteadyState,trend_coeff,aK] = DsgeSmoother(deep,gend,Y)
+  
+  stock_smooth(:,:,b) = alphahat;
+  if nvx
+    stock_innov(:,:,b)  = etahat;
+  end
+  if nvn
+    stock_error(:,:,b)  = epsilonhat;
+  end
+  if naK
+    stock_filter(:,:,:,b) = aK;
+  end
+
+  irun1 = irun1 + 1;
+  irun2 = irun2 + 1;
+  irun3 = irun3 + 1;
+  irun4 = irun4 + 1;
+
+  if irun1 == MAX_nsmoo | b == B
+    if b == B
+      stock_smooth = stock_smoo(:,:,1:irun1);
+    end
+    save([DirectoryName M_.fname '_smooth' int2str(ifil1)],'stock_smooth');
+    ifil1 = ifil1 + 1;
+    irun1 = 1;
+  end
+  
+  if nvx & (irun2 == MAX_inno | b == B)
+    if b == B
+      stock_innov = stock_innov(:,:,1:irun2);
+    end
+    save([DirectoryName M_.fname '_inno' int2str(ifil2)],'stock_inno');
+    ifil2 = ifil2 + 1;
+    irun2 = 1;
+  end
+    
+  if nvn & (irun3 == MAX_error | b == B)
+    if b == B
+      stock_error = stock_error(:,:,1:irun3);
+    end
+    save([DirectoryName M_.fname '_error' int2str(ifil3)],'stock_error');
+    ifil3 = ifil3 + 1;
+    irun3 = 1;
+  end
+    
+  if naK & (irun3 == MAX_naK | b == B)
+    if b == B
+      stock_filter = stock_filter(:,:,:,1:irun4);
+    end
+    save([DirectoryName M_.fname '_filter' int2str(ifil4)],'stock_filter');
+    ifil4 = ifil4 + 1;
+    irun4 = 1;
+  end
+    
+  waitbar(b/B,h);
+end
+close(h)
diff --git a/matlab/dynare_estimation.m b/matlab/dynare_estimation.m
index 3c6dd799cdcbbda008749bd73eb0c0f11fdb64e4..240152fcbc0f87c6bc6a07bfa7e4c87e5f9bc1ac 100644
--- a/matlab/dynare_estimation.m
+++ b/matlab/dynare_estimation.m
@@ -57,7 +57,10 @@ options_ = set_default_option(options_,'posterior_mode_estimation',1);
 options_ = set_default_option(options_,'MaxNumberOfBytes',1e6);
 options_ = set_default_option(options_,'xls_sheet','');
 options_ = set_default_option(options_,'xls_range','');
-options_ = set_default_option(options_.'filter_step_ahead',0);
+options_ = set_default_option(options_,'filter_step_ahead',0);
+if options_.filtered_vars ~= 0 & options_filter_step_ahead == 0
+  options_.filter_step_ahead = 1;
+end
 if options_.filter_step_ahead ~= 0
   options_.nk = max(options_.filter_step_ahead);
 else
diff --git a/matlab/dynare_m.exe b/matlab/dynare_m.exe
index 7f669242a868cd4e03c26ebb5a60520294423f1d..c56c9c38e85fb332850fc7db37eb6693ff177823 100755
Binary files a/matlab/dynare_m.exe and b/matlab/dynare_m.exe differ
diff --git a/matlab/metropolis_draw.m b/matlab/metropolis_draw.m
new file mode 100644
index 0000000000000000000000000000000000000000..2f864db8a2781444485a5dfb54a9406b7047a343
--- /dev/null
+++ b/matlab/metropolis_draw.m
@@ -0,0 +1,39 @@
+function xparam=metropolis_draw(init)
+  global options_ estim_params_
+  persistent mh_nblck NumberofDraws fname FirstLine FirstMhFile MAX_nruns
+  
+  if init
+    nvx  = estim_params_.nvx;
+    nvn  = estim_params_.nvn;
+    ncx  = estim_params_.ncx;
+    ncn  = estim_params_.ncn;
+    np   = estim_params_.np ;
+    npar = nvx+nvn+ncx+ncn+np;
+    DirectoryName = CheckPath('metropolis');
+    fname = [ DirectoryName '/' M_.fname]
+    load([ fname '_mh_history'])
+    FirstMhFile = record.KeepedDraws.FirstMhFile;
+    FirstLine = record.KeepedDraws.FirstLine; 
+    TotalNumberOfMhFiles = sum(record.MhDraws(:,2)); LastMhFile = TotalNumberOfMhFiles; 
+    TotalNumberOfMhDraws = sum(record.MhDraws(:,1));
+    NumberOfDraws = TotalNumberOfMhDraws-floor(options_.mh_drop* ...
+					       TotalNumberOfMhDraws);
+    MAX_nruns = ceil(options_.MaxNumberOfBytes/(npar+2)/8);
+    mh_nblck = options_.nblck;
+    return
+  end
+  
+  ChainNumber = ceil(rand*mh_nblck);
+  DrawNumber  = ceil(rand*NumberOfDraws);
+
+  if DrawNumber <= MAX_nruns-FirstLine+1
+    MhFilNumber = FirstMhFile;
+    MhLine = FirstLine+DrawNumber-1;
+  else
+    DrawNumber  = DrawNumber-(MAX_nruns-FirstLine+1);
+    MhFilNumber = FirstMhFile+ceil(DrawNumber/MAX_nruns); 
+    MhLine = DrawNumber-(MhFilNumber-FirstMhFile-1)*MAX_nruns;
+  end
+
+  load( [ fname '_mh' int2str(MhFilNumber) '_blck' int2str(ChainNumber) '.mat' ],'x2');
+  xparams = x2(MhLine,:);
diff --git a/matlab/set_all_parameters.m b/matlab/set_all_parameters.m
new file mode 100644
index 0000000000000000000000000000000000000000..61a7f9da52150f6817f8e167e4df3ecf47e18a45
--- /dev/null
+++ b/matlab/set_all_parameters.m
@@ -0,0 +1,56 @@
+function set_all_parameters(xparam1)
+  global estim_params_ M_
+  
+  nvx = estim_params_.nvx;
+  ncx = estim_params_.ncx;
+  nvn = estim_params_.nvn;
+  ncn = estim_params_.ncn;
+  np = estim_params_.np;
+  Sigma_e = M_.Sigma_e;
+  
+  if nvx
+    var_exo = estim_params_.var_exo;
+    for i=1:nvx
+      k =var_exo(i,1);
+      Sigma_e(k,k) = xparam1(i)^2;
+    end
+  end
+  
+  if ncx
+    offset = nvx+nvn;
+    corrx = estim_params_.corrx;
+    for i=1:ncx
+      k1 = corrx(i,1);
+      k2 = corrx(i,2);
+      Sigma_e(k1,k2) = xparam1(i+offset)*sqrt(Sigma_e_(k1,k1)*Sigma_e_(k2,k2));
+      Sigma_e(k2,k1) = Sigma_e_(k1,k2);
+    end
+  end
+  
+  if nvn
+    offset = nvx;
+    var_endo = estim_params_.var_endo;
+    for i=1:nvn
+      k = var_endo(i,1);
+      H(k,k) = xparam1(i+offset)^2;
+    end
+  end
+  
+  if ncx
+    offset = nvx+nvn+ncx;
+    corrn = estim_params_.corrn;
+    for i=1:ncx
+      k1 = corr(i,1);
+      k2 = corr(i,2);
+      H(k1,k2) = xparam1(i+offset)*sqrt(H(k1,k1)*H(k2,k2));
+      H(k2,k1) = H(k1,k2);
+    end
+  end
+  
+  if np
+    offset = offset+estim_params_.ncx+estim_params_.ncn;
+    M_.params(estim_params_.param_vals(:,1)) = deep(offset+1:end);
+  end
+  
+  M_.Sigma_e = Sigma_e;
+  M_.H = H;
\ No newline at end of file
diff --git a/matlab/set_parameters.m b/matlab/set_parameters.m
new file mode 100644
index 0000000000000000000000000000000000000000..62c30a851e80c21373a94beb457d2db746c3c518
--- /dev/null
+++ b/matlab/set_parameters.m
@@ -0,0 +1,33 @@
+function set_parameters(xparam1)
+  global estim_params_ M_
+  
+  nvx = estim_params_.nvx;
+  ncx = estim_params_.ncx;
+  np = estim_params_.np;
+  Sigma_e = M_.Sigma_e;
+
+  if nvx
+    var_exo = estim_params_.var_exo;
+    for i=1:nvx
+      k = var_exo(i,1);
+      Sigma_e(k,k) = xparam1(i)^2;
+    end
+  end
+  
+  if ncx
+    offset = nvx + nvn;
+    corrx = estim_params_.corrx;
+    for i=1:ncx
+      k1 = corrx(i,1);
+      k2 = corrx(i,2);
+      Sigma_e(k1,k2) = xparam1(i+offset)*sqrt(Sigma_e_(k1,k1)*Sigma_e_(k2,k2));
+      Sigma_e(k2,k1) = Sigma_e_(k1,k2);
+    end
+  end
+  
+  if np
+    offset = offset+estim_params_.ncx+estim_params_.ncn;
+    M_.params(estim_params_.param_vals(:,1)) = deep(offset+1:end);
+  end
+  
+  M_.Sigma_e = Sigma_e;