diff --git a/matlab/dynare_estimation.m b/matlab/dynare_estimation.m
index 8d839f33c..9b325bc10 100644
--- a/matlab/dynare_estimation.m
+++ b/matlab/dynare_estimation.m
@@ -1389,7 +1389,7 @@ if (~((any(bayestopt_.pshape > 0) & options_.mh_replic) | (any(bayestopt_.pshape
 end 
 
 if options_.forecast > 0 & options_.mh_replic == 0 & ~options_.load_mh_file 
-  forecast(var_list);
+  forecast(var_list_,'smoother');
 end
 
 pindx = estim_params_.param_vals(:,1);
diff --git a/matlab/forcst.m b/matlab/forcst.m
index 47dcefe88..d299da5ed 100644
--- a/matlab/forcst.m
+++ b/matlab/forcst.m
@@ -1,54 +1,77 @@
-function [yf,int_width]=forcst(dr,y0,k,var_list)
-  global M_  oo_ options_ 
-  
-  make_ex_;
-  yf = simult_(y0,dr,zeros(k,M_.exo_nbr),1);
-  nstatic = dr.nstatic;
-  npred = dr.npred;
-  nc = size(dr.ghx,2);
-  endo_nbr = M_.endo_nbr;
-  inv_order_var = dr.inv_order_var;
-  [A,B] = kalman_transition_matrix(dr,nstatic+(1:npred),1:nc,dr.transition_auxiliary_variables);
-  
-  nvar = size(var_list,1);
-  if nvar == 0
-    nvar = M_.endo_nbr;
-    ivar = [1:nvar];
-  else
-    ivar=zeros(nvar,1);
+function [yf,int_width]=forcst(dr,y0,horizon,var_list)
+
+% function [yf,int_width]=forecst(dr,y0,horizon,var_list)
+%   computes mean forecast for a given value of the parameters
+%   computes also confidence band for the forecast    
+%
+% INPUTS:
+%   dr:          structure containing decision rules
+%   y0:          initial values
+%   horizon:     nbr of periods to forecast
+%   var_list:    list of variables (character matrix)
+%
+% OUTPUTS:
+%   yf:          mean forecast
+%   int_width:   distance between upper bound and
+%                mean forecast
+%
+% SPECIAL REQUIREMENTS
+%    none
+%  
+% part of DYNARE, copyright Dynare Team (2003-2008)
+% Gnu Public License.
+    
+    global M_  oo_ options_ 
+    
+    make_ex_;
+    yf = simult_(y0,dr,zeros(horizon,M_.exo_nbr),1);
+    nstatic = dr.nstatic;
+    npred = dr.npred;
+    nc = size(dr.ghx,2);
+    endo_nbr = M_.endo_nbr;
+    inv_order_var = dr.inv_order_var;
+    [A,B] = kalman_transition_matrix(dr,nstatic+(1:npred),1:nc,dr.transition_auxiliary_variables);
+    
+    nvar = size(var_list,1);
+    if nvar == 0
+	nvar = M_.endo_nbr;
+	ivar = [1:nvar];
+    else
+	ivar=zeros(nvar,1);
+	for i=1:nvar
+	    i_tmp = strmatch(var_list(i,:),M_.endo_names,'exact');
+	    if isempty(i_tmp)
+		disp(var_list(i,:));
+		error (['One of the variable specified does not exist']) ;
+	    else
+		ivar(i) = i_tmp;
+	    end
+	end
+    end
+
+    ghx1 = dr.ghx(inv_order_var(ivar),:);
+    ghu1 = dr.ghu(inv_order_var(ivar),:);
+
+    sigma_u = B*M_.Sigma_e*B';
+    sigma_u1 = ghu1*M_.Sigma_e*ghu1';
+    sigma_y = 0;
+    
+    for i=1:horizon
+	sigma_y1 = ghx1*sigma_y*ghx1'+sigma_u1;
+	var_yf(i,:) = diag(sigma_y1)';
+	if i == horizon
+	    break
+	end
+	sigma_u = A*sigma_u*A';
+	sigma_y = sigma_y+sigma_u;
+    end
+
+    fact = qnorm((1-options_.conf_sig)/2,0,1);
+    
+    int_width = zeros(horizon,M_.endo_nbr);
     for i=1:nvar
-      i_tmp = strmatch(var_list(i,:),M_.endo_names,'exact');
-      if isempty(i_tmp)
-	disp(var_list(i,:));
-	error (['One of the variable specified does not exist']) ;
-      else
-	ivar(i) = i_tmp;
-      end
+	int_width(:,i) = fact*sqrt(var_yf(:,i));
     end
-  end
 
-  ghx1 = dr.ghx(inv_order_var(ivar),:);
-  ghu1 = dr.ghu(inv_order_var(ivar),:);
-
-  sigma_u = B*M_.Sigma_e*B';
-  sigma_u1 = ghu1*M_.Sigma_e*ghu1';
-  sigma_y = 0;
-  
-  for i=1:k
-    sigma_y1 = ghx1*sigma_y*ghx1'+sigma_u1;
-    var_yf(i,:) = diag(sigma_y1)';
-    if i == k
-      break
-    end
-    sigma_u = A*sigma_u*A';
-    sigma_y = sigma_y+sigma_u;
-  end
-
-  fact = qnorm((1-options_.conf_sig)/2,0,1);
-  
-  int_width = zeros(k,M_.endo_nbr);
-  for i=1:nvar
-    int_width(:,i) = fact*sqrt(var_yf(:,i));
-  end
-
-  yf = yf(ivar,:);
\ No newline at end of file
+    yf = yf(ivar,:);
+    
\ No newline at end of file
diff --git a/matlab/forecast.m b/matlab/forecast.m
index 844f55968..ab21b9301 100644
--- a/matlab/forecast.m
+++ b/matlab/forecast.m
@@ -1,49 +1,125 @@
-% Copyright (C) 2005 Michel Juillard
-%
-function forecast(var_list)
-  global options_ dr_ oo_ M_
-  
-  old_options = options_;
-  options_ = set_default_option(options_,'periods',40);
-  if options_.periods == 0
-    options_.periods = 40;
-  end
-  options_ = set_default_option(options_,'conf_sig',0.9);
-  
-  if size(oo_.endo_simul,2) < M_.maximum_lag
-    y0 = repmat(oo_.steady_state,M_.maximum_lag);
-  else
-    y0 = oo_.endo_simul(:,1:M_.maximum_lag);
-  end
-  
-  if M_.exo_det_nbr == 0
-    [yf,int_width] = forcst(oo_.dr,y0,options_.periods,var_list);
-  else
-    exo_det_length = size(oo_.exo_det_simul,1);
-    if options_.periods > exo_det_length
-      ex = zeros(options_.periods,M_.exo_nbr);
-      oo_.exo_det_simul = [ oo_.exo_det_simul;...
-		    repmat(oo_.exo_det_steady_state',...
-			   options_.periods- ... 
-			   exo_det_length,1)];
-			   %ex_det_length,1),1)];
-    elseif options_.periods < exo_det_length 
-      ex = zeros(exo_det_length,M_.exo_nbr); 
-    end
-    [yf,int_width] = simultxdet(y0,dr_,ex,oo_.exo_det_simul,...
-				options_.order,var_list);
-  end
-  
-  for i=1:M_.endo_nbr
-    eval(['oo_.forecast.Mean.' M_.endo_names(i,:) '= yf(' int2str(i) ',M_.maximum_lag+1:end)'';']);
-    eval(['oo_.forecast.HPDinf.' M_.endo_names(i,:) '= yf(' int2str(i) ',M_.maximum_lag+1:end)''-' ...
-		    ' int_width(:,' int2str(i) ');']);
-    eval(['oo_.forecast.HPDsup.' M_.endo_names(i,:) '= yf(' int2str(i) ',M_.maximum_lag+1:end)''+' ...
-		    ' int_width(:,' int2str(i) ');']);
-  end
+function forecast(var_list,task)
 
-  for i=1:M_.exo_det_nbr
-    eval(['oo_.forecast.Exogenous.' lgx_det_(i,:) '= M_.exo_det_simul(:,' int2str(i) ');']);
-  end
-  
-  options_ = old_options;
\ No newline at end of file
+% function forecast(var_list,task)
+%   computes mean forecast for a given value of the parameters
+%   computes also confidence band for the forecast    
+%
+% INPUTS
+%   var_list:    list of variables (character matrix)
+%   task:        indicates how to initialize the forecast
+%                either 'simul' or 'smoother'
+% OUTPUTS
+%   nothing is returned but the procedure saves output
+%   in oo_.forecast.Mean
+%      oo_.forecast.HPDinf
+%      oo_.forecast.HPDsup
+%
+% SPECIAL REQUIREMENTS
+%    none
+%  
+% part of DYNARE, copyright Dynare Team (2003-2008)
+% Gnu Public License.
+    
+    global options_ dr_ oo_ M_
+    
+    old_options = options_;
+
+    maximum_lag = M_.maximum_lag;
+    y_smoothed = oo_.SmoothedVariables;
+    options_ = set_default_option(options_,'periods',40);
+    if options_.periods == 0
+        options_.periods = 40;
+    end
+    horizon = options_.periods;
+    options_ = set_default_option(options_,'conf_sig',0.9);
+
+    endo_names = M_.endo_names;
+    if isempty(var_list)
+	varlist = endo_names;
+	i_var = 1:M_.endo_nbr;
+    else
+	i_var = [];
+	for i = 1:size(var_list)
+	    tmp = strmatch(var_list(i,:),endo_names,'exact');
+	    if isempty(tmp)
+		error([var_list(i,:) ' isn''t and endogenous variable'])
+	    end
+	    i_var = [i_var; tmp];
+	end
+    end
+    n_var = length(i_var);
+    
+    trend = 0;
+    switch task
+     case 'simul'
+      if size(oo_.endo_simul,2) < maximum_lag
+          y0 = repmat(oo_.steady_state,1,maximum_lag);
+      else
+          y0 = oo_.endo_simul(:,1:maximum_lag);
+      end
+     case 'smoother'
+      y0 = zeros(M_.endo_nbr,maximum_lag);
+      for i = 1:M_.endo_nbr
+          v_name = deblank(M_.endo_names(i,:));
+          y0(i,:) = y_smoothed.(v_name)(end-maximum_lag+1:end)+oo_.dr.ys(i);
+      end
+      gend = size(y0,2);
+      if isfield(oo_.Smoother,'TrendCoeffs')
+          var_obs = options_.varobs;
+          endo_names = M_.endo_names;
+          order_var = oo_.dr.order_var;
+          i_var_obs = [];
+          trend_coeffs = [];
+          for i=1:size(var_obs,1)
+              tmp = strmatch(var_obs(i,:),endo_names(i_var,:),'exact');
+              if ~isempty(tmp)
+                  i_var_obs = [ i_var_obs; tmp];
+                  trend_coeffs = [trend_coeffs; oo_.Smoother.TrendCoeffs(i)];
+              end
+          end         
+          trend = trend_coeffs*(gend+(1-M_.maximum_lag:horizon));
+      end
+     otherwise
+      error('Wrong flag value')
+    end 
+    
+    if M_.exo_det_nbr == 0
+        [yf,int_width] = forcst(oo_.dr,y0,options_.periods,var_list);
+    else
+        exo_det_length = size(oo_.exo_det_simul,1);
+        if options_.periods > exo_det_length
+            ex = zeros(options_.periods,M_.exo_nbr);
+            oo_.exo_det_simul = [ oo_.exo_det_simul;...
+                    repmat(oo_.exo_det_steady_state',...
+                           options_.periods- ... 
+                           exo_det_length,1)];
+            %ex_det_length,1),1)];
+        elseif options_.periods < exo_det_length 
+            ex = zeros(exo_det_length,M_.exo_nbr); 
+        end
+        [yf,int_width] = simultxdet(y0,dr_,ex,oo_.exo_det_simul,...
+                                    options_.order,var_list);
+    end
+
+    if ~isscalar(trend)
+        yf(i_var_obs,:) = yf(i_var_obs,:) + trend;
+    end
+    
+    for i=1:n_var
+        eval(['oo_.forecast.Mean.' var_list(i,:) '= yf(' int2str(i) ',maximum_lag+1:end)'';']);
+        eval(['oo_.forecast.HPDinf.' var_list(i,:) '= yf(' int2str(i) ',maximum_lag+1:end)''-' ...
+              ' int_width(:,' int2str(i) ');']);
+        eval(['oo_.forecast.HPDsup.' var_list(i,:) '= yf(' int2str(i) ',maximum_lag+1:end)''+' ...
+              ' int_width(:,' int2str(i) ');']);
+    end
+
+    for i=1:M_.exo_det_nbr
+        eval(['oo_.forecast.Exogenous.' lgx_det_(i,:) '= M_.exo_det_simul(:,' int2str(i) ');']);
+    end
+    
+    if options_.nograph == 0
+        forecast_graphs(var_list);
+    end
+    
+    options_ = old_options;
+    
diff --git a/matlab/forecast_graphs.m b/matlab/forecast_graphs.m
new file mode 100644
index 000000000..99c0eec06
--- /dev/null
+++ b/matlab/forecast_graphs.m
@@ -0,0 +1,77 @@
+function forecast_graphs(var_list)
+    
+    global M_ oo_ options_
+    
+    nc = 4;
+    nr = 3;
+    exo_nbr = M_.exo_nbr;
+    endo_names = M_.endo_names;
+    fname = M_.fname;
+% $$$     varobs = options_.varobs;
+% $$$     y = oo_.SmoothedVariables;
+% $$$     ys = oo_.dr.ys;
+% $$$     gend = size(y,2);
+    yf = oo_.forecast.Mean;
+    hpdinf = oo_.forecast.HPDinf;
+    hpdsup = oo_.forecast.HPDsup;
+    if isempty(var_list)
+	varlist = endo_names;
+	i_var = 1:M_.endo_nbr;
+    else
+	i_var = [];
+	for i = 1:size(var_list)
+	    tmp = strmatch(var_list(i,:),endo_names,'exact');
+	    if isempty(tmp)
+		error([var_list(i,:) ' isn''t and endogenous variable'])
+	    end
+	    i_var = [i_var; tmp];
+	end
+    end
+    nvar = length(i_var);
+    
+% $$$     % build trend for smoothed variables if necessary
+% $$$     trend = zeros(size(varobs,1),10);
+% $$$     if isfield(oo_.Smoother,'TrendCoeffs')
+% $$$         trend_coeffs = oo_.Smoother.TrendCoeffs;
+% $$$         trend = trend_coeffs*(gend-9:gend);
+% $$$     end
+
+    % create subdirectory <fname>/graphs if id doesn't exist
+    if ~exist([fname '/graphs'])
+	mkdir(fname,'graphs');
+    end
+    
+    m = 1;
+    n_fig = 1;
+    figure('Name','Forecasts (I)')
+    for j= 1:nvar
+	if m > nc*nr; 
+	    eval(['print -depsc ' fname '/graphs/forcst' int2str(n_fig) ...
+		 '.eps;'])
+	    n_fig =n_fig+1;
+	    eval(['figure(''Name'',''Forecast (' int2str(n_fig) ')'');']);
+	    m = 1;
+	end
+	subplot(nr,nc,m);
+	vn = deblank(endo_names(i_var(j),:));
+	obs = 0;
+% $$$         k = strmatch(vn,varobs,'exact'); 
+% $$$ 	if ~isempty(k)
+% $$$ 	    yy = y.(vn)(end-9:end) + repmat(ys(i_var(j)),10,1)+trend(k,:)';
+% $$$ 	    plot(yy);
+% $$$ 	    hold on
+% $$$ 	    obs = 10;
+% $$$ 	end
+	plot([NaN(obs,1); yf.(vn)]);
+	hold on
+	plot([NaN(obs,1); hpdinf.(vn)]);
+	hold on
+	plot([NaN(obs,1); hpdsup.(vn)]);
+	title(vn,'Interpreter','none');
+	hold off
+	m = m + 1;
+    end
+    
+    if m > 1
+	eval(['print -deps ' fname '/graphs/forcst' int2str(n_fig) '.eps;'])
+    end
\ No newline at end of file