diff --git a/matlab/dynare_estimation.m b/matlab/dynare_estimation.m
index 44cae711dce11966bed303fdbce8d2db281004ad..fd6098f679b330f26ec52ae8329e661b01df287a 100644
--- a/matlab/dynare_estimation.m
+++ b/matlab/dynare_estimation.m
@@ -72,25 +72,17 @@ n_varobs = size(options_.varobs,1);
 dr = set_state_space([]);
 
 % load data
-if exist(options_.datafile)
-  instr = options_.datafile;
-else
-  instr = ['load ' options_.datafile];
-end
-
-eval(instr);
+rawdata = read_variables(options_.datafile,options_.varobs,[]);
 
-rawdata = [];
 k = [];
 k1 = [];
 for i=1:n_varobs
-  rawdata = [rawdata eval(deblank(options_.varobs(i,:)))];
   k = [k strmatch(deblank(options_.varobs(i,:)),lgy_(dr.order_var,:), ...
           'exact')];
   k1 = [k1 strmatch(deblank(options_.varobs(i,:)),lgy_, 'exact')];
 end
 
-% check initial parameter values
+% check initial parameter values and creates bayestopt_
 pnames=['     ';'beta ';'gamm ';'norm ';'invg ';'unif ';'invg2'];
 [xparam1,estim_params_,bayestopt_,lb,ub]=set_prior(estim_params_);
 if any(bayestopt_.pshape > 0)
diff --git a/matlab/dynare_m.exe b/matlab/dynare_m.exe
index 824ac796bdd9a8e03ef70ef211f8c9d24bfb30d1..ba9cb681080d8540772ff2a36da25d2c39bc8703 100644
Binary files a/matlab/dynare_m.exe and b/matlab/dynare_m.exe differ
diff --git a/matlab/forcst.m b/matlab/forcst.m
index 86bc921bc3edde5bfe3aa934c5d1f0e93487ff94..afabf587287890f794b63a0e60b59c20998f679d 100644
--- a/matlab/forcst.m
+++ b/matlab/forcst.m
@@ -1,4 +1,4 @@
-function [yf,var_yf]=forcst(dr,y0,k,var_list)
+function [yf,var_yf] = forcst(dr,y0,k,var_list)
   global endo_nbr exo_nbr ykmin_ Sigma_e_ ex_ options_ lgy_
   
   options_.periods = k;
diff --git a/matlab/forecast.m b/matlab/forecast.m
new file mode 100644
index 0000000000000000000000000000000000000000..ec4ae16a5b36a7e28ae0c710151cee8af5e047ef
--- /dev/null
+++ b/matlab/forecast.m
@@ -0,0 +1,39 @@
+% Copyright (C) 2005 Michel Juillard
+%
+function forecast(var_list)
+  global options_ dr_ ys_ endo_nbr exo_nbr exo_det_nbr ykmin_ y_ ex_det_ 
+  global lgy_ lgx_det_ oo_ exe_det_
+  
+  old_options = options_;
+  options_ = set_default_option(options_,'periods',40);
+  if options_.periods == 0
+    options_.periods = 40;
+  end
+  
+  if size(y_,2) < ykmin_
+    y0 = repmat(ys_,ykmin_);
+  else
+    y0 = y_(:,1:ykmin_);
+  end
+  
+  if exo_det_nbr == 0
+    [yf,var_yf] = forcst(dr_,y0,options_.periods,var_list);
+  else
+    ex = zeros(options_.periods,exo_nbr);
+    if options_.periods > size(ex_det_,1)
+      ex_det_ = [ ex_det_; repmat(exe_det_',options_.periods- ...
+				  size(ex_det_,1),1)];
+    end
+    [yf,var_yf] = simultxdet(y0,dr_,ex,ex_det_,options_.order, ...
+				  var_list);
+  end
+  
+  for i=1:endo_nbr
+    eval(['oo_.forecast.Mean.' lgy_(i,:) '= yf(' int2str(i) ',:)'';']);
+  end
+
+  for i=1:exo_det_nbr
+    eval(['oo_.forecast.Mean.' lgx_det_(i,:) '= ex_det_(:,' int2str(i) ');']);
+  end
+  
+  options_ = old_options;
\ No newline at end of file
diff --git a/matlab/read_variables.m b/matlab/read_variables.m
new file mode 100644
index 0000000000000000000000000000000000000000..38d821291d59b08ea4019c5aa581071b79ac144a
--- /dev/null
+++ b/matlab/read_variables.m
@@ -0,0 +1,25 @@
+% Copyright (C) 2005 Michel Juillard
+%
+% all local variables have complicated names in order to avoid name
+% conflicts with possible user variable names
+
+function dyn_data_01=read_variables(file_name_01,var_names_01,dyn_data_01)
+  
+  dyn_size_01 = size(dyn_data_01,1);
+  
+  if exist(file_name_01)
+    dyn_instr_01 = file_name_01;
+  else
+    dyn_instr_01 = ['load ' file_name_01];
+  end
+  
+  eval(dyn_instr_01);
+
+  
+  for dyn_i_01=1:size(var_names_01,1)
+    dyn_tmp_01 = eval(var_names_01(dyn_i_01,:));
+    if length(dyn_tmp_01) > dyn_size_01 & dyn_size_01 > 0
+      error('data size is too large')
+    end
+    dyn_data_01(:,dyn_i_01) = dyn_tmp_01;
+  end
\ No newline at end of file
diff --git a/matlab/shocks_file.m b/matlab/shocks_file.m
new file mode 100644
index 0000000000000000000000000000000000000000..76a73a61b179926b9810de0a06643f04017731b4
--- /dev/null
+++ b/matlab/shocks_file.m
@@ -0,0 +1,44 @@
+% Copyright (C) 2001 Michel Juillard
+%
+function []=shocks_file(ms_flag)
+  global options_ M_ exo_nbr exo_det_nbr lgx_ lgx_det_ ex_ exe_ ex_det_ exe_det_
+  
+  if isfield('options_','shocks_file')
+    if exo_det_nbr > 0
+      x = read_variables(options_.shocks_file,lgx_det_);
+      n = size(x,1);
+      M_.ex_det_length = n;
+      if size(ex_det_,1) >= n
+	if ms_flag == 1
+	  ex_det_(1:n,:) = ex_det_(1:n,:).*x;
+	else
+	  ex_det_(1:n,:) = x;
+	end
+      else
+	if ms_flag == 1
+	  ex_det_ = ones(n,1)*exe_det_';
+	  ex_det_ = ex_det_.*x;
+	else
+	  ex_det_ = x;
+	end
+      end
+    else
+      x = read_variables(options_.shocks_file,lgx_det_);
+      n = size(x,1);
+      if size(ex_,1) >= n
+	if ms_flag == 1
+	  ex_(1:n,:) = ex_(1:n,:).*x;
+	else
+	  ex_(1:n,:) = x;
+	end
+      else
+	if ms_flag == 1
+	  ex_ = ones(n,1)*exe_';
+	  ex_ = ex_.*x;
+	else
+	  ex_ = x;
+	end
+      end
+    end
+  end
+  
\ No newline at end of file
diff --git a/matlab/simultxdet.m b/matlab/simultxdet.m
new file mode 100644
index 0000000000000000000000000000000000000000..0bd531e9757c267dda63239d40c003bc0e8d1a21
--- /dev/null
+++ b/matlab/simultxdet.m
@@ -0,0 +1,136 @@
+% Copyright (C) 2005 Michel Juillard
+%
+
+function [y_,var_yf]=simultexdet(y0,dr,ex_,ex_det, iorder,var_list)
+  global endo_nbr ykmin_ xkmin_ it_ options_ iy_ M_ exe_det_ Sigma_e_ lgy_
+
+  iter = size(ex_,1)-xkmin_;
+  nx = size(dr.ghu,2);
+  y_ = zeros(size(y0,1),iter+ykmin_);
+  y_(:,1:ykmin_) = y0;
+  k1 = [ykmin_:-1:1];
+  k2 = dr.kstate(find(dr.kstate(:,2) <= ykmin_+1),[1 2]);
+  k2 = k2(:,1)+(ykmin_+1-k2(:,2))*endo_nbr;
+  k3 = iy_(1:ykmin_,:)';
+  k3 = find(k3(:));
+  k4 = dr.kstate(find(dr.kstate(:,2) < ykmin_+1),[1 2]);
+  k4 = k4(:,1)+(ykmin_+1-k4(:,2))*endo_nbr;
+  
+  old_options = options_;
+  options_ = set_default_option(options_,'simul_algo',0);
+  options_ = set_default_option(options_,'conf_sig',0.9);
+  
+  if options_.simul_algo == 1
+    o1 = dr.nstatic+1;
+    o2 = dr.nstatic+dr.npred;
+    o3 = o2-dr.nboth+1;
+    [junk, k5] = sort(dr.order_var(o1:o2));
+    [junk, k6] = sort(dr.order_var(o3:end));
+  end
+  
+  nvar = size(var_list,1);
+  if nvar == 0
+    nvar = endo_nbr;
+    ivar = [1:nvar];
+  else
+    ivar=zeros(nvar,1);
+    for i=1:nvar
+      i_tmp = strmatch(var_list(i,:),lgy_,'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
+
+  if iorder == 1
+    for i = ykmin_+1: iter+ykmin_
+      tempx1 = y_(dr.order_var,k1);
+      tempx2 = tempx1-repmat(dr.ys(dr.order_var),1,ykmin_);
+      tempx = tempx2(k2);
+      if options_.simul_algo == 0
+	y_(dr.order_var,i) = dr.ys(dr.order_var)+dr.ghx*tempx+dr.ghu* ...
+	    ex_(i+xkmin_-ykmin_,:)';
+	for j=1:min(iter+ykmin_-i-M_.ex_det_length+1,M_.ex_det_length)
+	  y_(dr.order_var,i) = y_(dr.order_var,i) + dr.ghud{j}*(ex_det(i+j-1,:)'-exe_det_');
+	end
+      elseif options_.simul_algo == 1
+	it_ = i;
+	m = dr.ys(dr.order_var);
+	[y_(:,i), check] = dynare_solve('ff_simul1',y_(:,i-1),tempx1(k3), ...
+					m(o3:end),tempx(k4),o1,o2,o3,k6);
+      end
+	
+      k1 = k1+1;
+    end
+  elseif iorder == 2
+    for i = ykmin_+1: iter+ykmin_
+      tempx1 = y_(dr.order_var,k1);
+      tempx2 = tempx1-repmat(dr.ys(dr.order_var),1,ykmin_);
+      tempx = tempx2(k2);
+      tempu = ex_(i+xkmin_-ykmin_,:)';
+      tempuu = kron(tempu,tempu);
+      if options_.simul_algo == 0
+	tempxx = kron(tempx,tempx);
+	tempxu = kron(tempx,tempu);
+	y_(dr.order_var,i) = dr.ys(dr.order_var)+dr.ghs2/2+dr.ghx*tempx+ ...
+	    dr.ghu*tempu+0.5*(dr.ghxx*tempxx+dr.ghuu*tempuu)+dr.ghxu* ...
+	    tempxu;
+	for j=1:min(iter+ykmin_-i-M_.ex_det_length+1,M_.ex_det_length)
+	  tempud = ex_det(i+j-1,:)'-exe_det_';
+	  tempudud = kron(tempud,tempud);
+	  tempxud = kron(tempx,tempud);
+	  tempuud = kron(tempu,tempud);
+	  y_(dr.order_var,i) = y_(dr.order_var,i) + dr.ghud{j}*tempud + ...
+	      dr.ghxud{j}*tempxud + dr.ghuud{j}*tempuud + ...
+	      0.5*dr.ghudud{j,j}*tempudud;
+	  for k=1:j-1
+	    tempudk = ex_det(i+k-1,:)'-exe_det_';
+	    tempududk = kron(tempudk,tempud);
+	    y_(dr.order_var,i) = y_(dr.order_var,i) + ...
+		dr.ghudud{k,j}*tempududk;
+	  end
+	end
+      elseif options_.simul_algo == 1
+	it_ = i;
+	m = dr.ys(dr.order_var)+dr.ghs2/2;
+	tempx1 = y_(:,k1);
+	[y_(:,i), check] = dynare_solve('ff_simul2',y_(:,i-1),tempx1(k3), ...
+					m(o3:end),tempx(k4),o1,o2,o3,k6);
+      end
+      k1 = k1+1;
+    end
+  end
+
+  [A,B] = kalman_transition_matrix(dr);
+  
+  sigma_u = B*Sigma_e_*B';
+  sigma_y = 0;
+  
+  for i=1:iter
+    sigma_y = sigma_y+sigma_u;
+    var_yf(i,dr.order_var) = diag(sigma_y(1:endo_nbr,1:endo_nbr))';
+    if i == iter
+      break
+    end
+    sigma_u = A*sigma_u*A';
+  end
+
+  fact = qnorm((1-options_.conf_sig)/2,0,1);
+  
+  for i=1:nvar
+    my_subplot(i,nvar,2,3,'Forecasts');
+    
+    interval_width = fact*sqrt(var_yf(:,ivar(i)));
+    
+    plot([-ykmin_+1:0],y0(ivar(i),1:ykmin_),'b-',...
+	 [1:iter],y_(ivar(i),ykmin_+1:end),'g-',...
+	 [1:iter],y_(ivar(i),ykmin_+1:end)'+interval_width,'g:',...
+	 [1:iter],y_(ivar(i),ykmin_+1:end)'-interval_width,'g:',...
+	 [1 iter],repmat(dr.ys(ivar(i)),1,2),'r-');
+    title(lgy_(ivar(i),:));
+  end
+
+  options_ = old_options;
\ No newline at end of file