diff --git a/matlab/forcst.m b/matlab/forcst.m
index 4e4a80037..47dcefe88 100644
--- a/matlab/forcst.m
+++ b/matlab/forcst.m
@@ -1,74 +1,54 @@
-function [yf,var_yf]=forcst(dr,y0,k,m)
-  global M_  oo_ options_
+function [yf,int_width]=forcst(dr,y0,k,var_list)
+  global M_  oo_ options_ 
   
-  options_.periods = k;
   make_ex_;
-  yf = simult_(y0,dr,oo_.exo_simul(1:k,:),1);
-
+  yf = simult_(y0,dr,zeros(k,M_.exo_nbr),1);
   nstatic = dr.nstatic;
   npred = dr.npred;
-  %j = find(dr.kstate(dr.kae,2) <= ykmin_+1);
-  j = find(dr.kstate(dr.kae,2) <= M_.maximum_lag+1);
-  kae = dr.kae(j);
-  nh = size(dr.ghx,2);
-  hx = dr.ghx(nstatic+1:nstatic+npred,:);
-  hu = dr.ghu(nstatic+1:nstatic+npred,:);
-  if ~isempty(kae)
-    n = length(kae);
-    tmp = sparse([1:n]',kae-sum(dr.kstate(:,2)>M_.maximum_lag+1),ones(n,1),n,nh);
-    hx = [hx; tmp];
-    hu = [hu; zeros(n,M_.exo_nbr)];
-  end
-  gx  = [];
-  k2  = [];
-  if nstatic > 0
-    gx = dr.ghx(1:nstatic,:);
-    k2 = [1:nstatic]';
-  end
-  if size(dr.ghx,1) > nstatic+npred
-    gx = [gx; dr.ghx(nstatic+npred+1:end,:)];
-    k2 = [k2; [nstatic+npred+1:size(dr.ghx,1)]'];
-  end
+  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);
   
-  k1 = dr.order_var([nstatic+1:nstatic+npred]);
-  k2 = dr.order_var(k2);
-  
-  sigma_u = hu*M_.Sigma_e*hu';
-  sigma_y1 = 0;
-  var_yf = zeros(k,M_.endo_nbr);
-  
-
-  if isempty(k2)
-    for i=1:k
-      sigma_y1 = sigma_y1+sigma_u;
-      var_yf(i,k1) = diag(sigma_y1(1:npred,1:npred))';
-      if i == k
-	break
-      end
-      sigma_u = hx*sigma_u*hx';
-    end
+  nvar = size(var_list,1);
+  if nvar == 0
+    nvar = M_.endo_nbr;
+    ivar = [1:nvar];
   else
-    for i=1:k
-      sigma_y1 = sigma_y1+sigma_u;
-      var_yf(i,k1) = diag(sigma_y1(1:npred,1:npred))';
-      sigma_y2 = gx*sigma_y1*gx';
-      var_yf(i,k2) = diag(sigma_y2)';
-      if i == k
-	  break
+    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
-      sigma_u = hx*sigma_u*hx';
     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
diff --git a/matlab/forcst_unc.m b/matlab/forcst_unc.m
index 74c2141a6..e01d19fd0 100644
--- a/matlab/forcst_unc.m
+++ b/matlab/forcst_unc.m
@@ -23,6 +23,10 @@ function forcst_unc(y0,var_list)
   % setting up estim_params_
   [xparam1,estim_params_,bayestopt_,lb,ub] = set_prior(estim_params_);
 
+  options_.TeX = 0;
+  options_.nograph = 0;
+  plot_priors;
+  
   % workspace initialization
   if isempty(var_list)
     var_list = M_.endo_names;
@@ -34,6 +38,7 @@ function forcst_unc(y0,var_list)
   exo_nbr = M_.exo_nbr;
   replic = options_.replic;
   order = options_.order;
+  maximum_lag = M_.maximum_lag;
   %  params = prior_draw(1);
   params = rndprior(bayestopt_);
   set_parameters(params);
@@ -103,15 +108,29 @@ function forcst_unc(y0,var_list)
     k2(2) = 1;
   end
 
+  % compute shock uncertainty around forecast with mean prior
+  set_parameters(bayestopt_.pmean);
+  [dr,info] = resol(oo_.steady_state,0);
+  [yf3,yf3_intv] = forcst(dr,y0,periods,var_list);
+  yf3_1 = yf3'-[zeros(maximum_lag,n); yf3_intv];
+  yf3_2 = yf3'+[zeros(maximum_lag,n); yf3_intv];
+  
   % graphs
   
-  dynare_graph_init('Forecasts',n,{'b-' 'g-' 'g-' 'r-' 'r-'});
+  dynare_graph_init('Forecasts type I',n,{'b-' 'g-' 'g-' 'r-' 'r-'});
   for i=1:n
     dynare_graph([yf_mean(:,i) squeeze(yf1(:,i,k1)) squeeze(yf2(:,i,k2))],...
 		 var_list(i,:));
   end
   dynare_graph_close;
   
+  dynare_graph_init('Forecasts type II',n,{'b-' 'k-' 'k-' 'r-' 'r-'});
+  for i=1:n
+    dynare_graph([yf_mean(:,i) yf3_1(:,i) yf3_2(:,i) squeeze(yf2(:,i,k2))],...
+		 var_list(i,:));
+  end
+  dynare_graph_close;
+  
   
   % saving results
   save_results(yf_mean,'oo_.forecast.mean.',var_list);
diff --git a/matlab/simult.m b/matlab/simult.m
index 3e26466c0..3f09bc8d9 100644
--- a/matlab/simult.m
+++ b/matlab/simult.m
@@ -7,6 +7,9 @@ global  it_ means_
 
 order = options_.order;
 replic = options_.replic;
+if replic == 0
+  replic = 1;
+end
 seed = options_.simul_seed;
 options_.periods = options_.periods;