v4: added forecast, updated dynare_m.exe

git-svn-id: https://www.dynare.org/svn/dynare/dynare_v4@1812 ac1d8469-bf42-47a9-8791-bf33cf982152

matlab/PosteriorFilterSmootherAndForecast.m

@@ -94,12 +94,12 @@ irun5 = 1;
 irun6 = 1;
 irun7 = 1;
 ifil1 = 0;
-ifil2 = 1;
-ifil3 = 1;
-ifil4 = 1;
-ifil5 = 1;
-ifil6 = 1;
-ifil7 = 1;
+ifil2 = 0;
+ifil3 = 0;
+ifil4 = 0;
+ifil5 = 0;
+ifil6 = 0;
+ifil7 = 0;
 
 h = waitbar(0,'Bayesian smoother...');
 
@@ -110,14 +110,15 @@ if options_.smoother
   stock_smooth = zeros(endo_nbr,gend,MAX_nsmoo);
   stock_innov  = zeros(exo_nbr,gend,B);
   stock_error = zeros(nvobs,gend,MAX_nerro);
-  if options_.filter_step_ahead
+end
+if options_.filter_step_ahead
     stock_filter = zeros(naK,endo_nbr,gend+options_.filter_step_ahead(end),MAX_naK);
-  end
-  if options_.forecast
-    stock_forcst_mean = zeros(endo_nbr,horizon+maxlag,MAX_forcst1);
-    stock_forcst_total = zeros(endo_nbr,horizon+maxlag,MAX_forcst2);
-  end
 end
+if options_.forecast
+    stock_forcst_mean = zeros(endo_nbr,horizon+maxlag,MAX_nforc1);
+    stock_forcst_total = zeros(endo_nbr,horizon+maxlag,MAX_nforc2);
+end
+
 for b=1:B
   %deep = GetOneDraw(NumberOfDraws,FirstMhFile,LastMhFile,FirstLine,MAX_nruns,DirectoryName);
   [deep, logpo] = GetOneDraw(type);
@@ -156,7 +157,7 @@ for b=1:B
     yf(:,IdObs) = yf(:,IdObs)+(gend+[1-maxlag:horizon]')*trend_coeff';
     if options_.loglinear == 1
       yf = yf+repmat(log(SteadyState'),horizon+maxlag,1);
-      yf = exp(yf);
+%      yf = exp(yf);
     else
       yf = yf+repmat(SteadyState',horizon+maxlag,1);
     end
@@ -169,13 +170,13 @@ for b=1:B
 					   trend_coeff',[1,1,1]);
     if options_.loglinear == 1
       yf1 = yf1 + repmat(log(SteadyState'),[horizon+maxlag,1,1]);
-      yf1 = exp(yf1);
+%      yf1 = exp(yf1);
     else
       yf1 = yf1 + repmat(SteadyState',[horizon+maxlag,1,1]);
     end
 
-    stock_forcst_mean(:,:,irun6) = yf;
-    stock_forcst_total(:,:,irun7) = yf1;
+    stock_forcst_mean(:,:,irun6) = yf';
+    stock_forcst_total(:,:,irun7) = yf1';
   end
   
   irun1 = irun1 + 1;
@@ -195,44 +196,44 @@ for b=1:B
   
   if nvx & (irun2 > MAX_ninno | b == B)
     stock = stock_innov(:,:,1:irun2-1);
-    save([DirectoryName '/' M_.fname '_inno' int2str(ifil2)],'stock');
     ifil2 = ifil2 + 1;
+    save([DirectoryName '/' M_.fname '_inno' int2str(ifil2)],'stock');
     irun2 = 1;
   end
     
   if nvn & (irun3 > MAX_error | b == B)
     stock = stock_error(:,:,1:irun3-1);
-    save([DirectoryName '/' M_.fname '_error' int2str(ifil3)],'stock');
     ifil3 = ifil3 + 1;
+    save([DirectoryName '/' M_.fname '_error' int2str(ifil3)],'stock');
     irun3 = 1;
   end
     
   if naK & (irun4 > MAX_naK | b == B)
     stock = stock_filter(:,:,:,1:irun4-1);
-    save([DirectoryName '/' M_.fname '_filter' int2str(ifil4)],'stock');
     ifil4 = ifil4 + 1;
+    save([DirectoryName '/' M_.fname '_filter' int2str(ifil4)],'stock');
     irun4 = 1;
   end
     
   if irun5 > MAX_nruns | b == B
     stock = stock_param(1:irun5-1,:);
-    save([DirectoryName '/' M_.fname '_param' int2str(ifil5)],'stock','stock_logpo','stock_ys');
     ifil5 = ifil5 + 1;
+    save([DirectoryName '/' M_.fname '_param' int2str(ifil5)],'stock','stock_logpo','stock_ys');
     irun5 = 1;
   end
 
   if horizon & (irun6 > MAX_nforc1 | b == B)
     stock = stock_forcst_mean(:,:,1:irun6-1);
-    save([DirectoryName '/' M_.fname '_forc_mean' int2str(ifil6)],'stock');
     ifil6 = ifil6 + 1;
+    save([DirectoryName '/' M_.fname '_forc_mean' int2str(ifil6)],'stock');
     irun6 = 1;
   end
 
   if horizon & (irun7 > MAX_nforc2 |  b == B)
     stock = stock_forcst_total(:,:,1:irun7-1);
+    ifil7 = ifil7 + 1;
     save([DirectoryName '/' M_.fname '_forc_total' int2str(ifil7)],'stock');
-    ifil6 = ifil6 + 1;
-    irun6 = 1;
+    irun7 = 1;
   end
 
   waitbar(b/B,h);
@@ -243,6 +244,17 @@ stock_gend=gend;
 stock_data=Y;
 save([DirectoryName '/' M_.fname '_data'],'stock_gend','stock_data');
 
-pm3(endo_nbr,gend,ifil1,B,'Smoothed variables',...
-    M_.endo_names(SelecVariables),M_.endo_names,'tit_tex',M_.endo_names,...
-    'names2','smooth',[M_.fname '/metropolis'],'_smooth')
\ No newline at end of file
+if options_.smoother
+    pm3(endo_nbr,gend,ifil1,B,'Smoothed variables',...
+	M_.endo_names(SelecVariables),M_.endo_names,'tit_tex',M_.endo_names,...
+	'names2','smooth',[M_.fname '/metropolis'],'_smooth')
+end
+
+if options_.forecast
+    pm3(endo_nbr,horizon+maxlag,ifil6,B,'Forecasted variables (mean)',...
+	M_.endo_names(SelecVariables),M_.endo_names,'tit_tex',M_.endo_names,...
+	'names2','smooth',[M_.fname '/metropolis'],'_forc_mean')
+    pm3(endo_nbr,horizon+maxlag,ifil6,B,'Forecasted variables (total)',...
+	M_.endo_names(SelecVariables),M_.endo_names,'tit_tex',M_.endo_names,...
+	'names2','smooth',[M_.fname '/metropolis'],'_forc_total')
+end

matlab/forcst2a.m

+function yf=forcst2a(y0,dr,e)
+  global M_ options_
+  
+  Sigma_e_ = M_.Sigma_e;
+  endo_nbr = M_.endo_nbr;
+  exo_nbr = M_.exo_nbr;
+  ykmin_ = M_.maximum_endo_lag;
+  
+  horizon = size(e,1);
+  options_ = set_default_option(options_,'simul_seed',0);
+  order = options_.order;
+
+  k1 = [ykmin_:-1:1];
+  k2 = dr.kstate(find(dr.kstate(:,2) <= ykmin_+1),[1 2]);
+  k2 = k2(:,1)+(ykmin_+1-k2(:,2))*endo_nbr;
+
+  yf = zeros(horizon+ykmin_,endo_nbr);
+  yf(1:ykmin_,:) = y0';
+  
+  j = ykmin_*endo_nbr;
+  for i=ykmin_+(1:horizon)
+    tempx = yf(k1,:)';
+    yf(i,:) = tempx(k2)'*dr.ghx';
+    k1 = k1+1;
+  end
+  
+   yf(:,dr.order_var) = yf;
+   

matlab/forecast.m

-function forecast(var_list,task)
+function info = forecast(var_list,task)
 
 % function forecast(var_list,task)
 %   computes mean forecast for a given value of the parameters
@@ -22,10 +22,11 @@ function forecast(var_list,task)
     
     global options_ dr_ oo_ M_
     
+    info = 0;
+    
     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;
@@ -58,6 +59,7 @@ function forecast(var_list,task)
           y0 = oo_.endo_simul(:,1:maximum_lag);
       end
      case 'smoother'
+      y_smoothed = oo_.SmoothedVariables;
       y0 = zeros(M_.endo_nbr,maximum_lag);
       for i = 1:M_.endo_nbr
           v_name = deblank(M_.endo_names(i,:));

preprocessor/ComputingTasks.cc

@@ -164,6 +164,32 @@ StochSimulStatement::writeOutput(ostream &output, const string &basename) const
   output << "info = stoch_simul(var_list_);\n";
 }
 
+ForecastStatement::ForecastStatement(const SymbolList &symbol_list_arg,
+                                         const OptionsList &options_list_arg) :
+  symbol_list(symbol_list_arg),
+  options_list(options_list_arg)
+{
+}
+
+void
+ForecastStatement::checkPass(ModFileStructure &mod_file_struct)
+{
+  mod_file_struct.stoch_simul_or_similar_present = true;
+
+  // Fill in option_order of mod_file_struct
+  OptionsList::num_options_type::const_iterator it = options_list.num_options.find("order");
+  if (it != options_list.num_options.end())
+    mod_file_struct.order_option = max(mod_file_struct.order_option,atoi(it->second.c_str()));
+}
+
+void
+ForecastStatement::writeOutput(ostream &output, const string &basename) const
+{
+  options_list.writeOutput(output);
+  symbol_list.writeOutput("var_list_", output);
+  output << "info = forecast(var_list_,'simul');\n";
+}
+
 RamseyPolicyStatement::RamseyPolicyStatement(const SymbolList &symbol_list_arg,
                                              const OptionsList &options_list_arg) :
   symbol_list(symbol_list_arg),

preprocessor/DynareBison.yy

@@ -189,6 +189,7 @@ statement : declaration
           | bvar_forecast
           | dynare_sensitivity
           | homotopy_setup
+          | forecast
           ;
 
 declaration : parameters
@@ -1276,6 +1277,20 @@ homotopy_item : NAME COMMA expression COMMA expression ';'
                 { driver.homotopy_val($1, NULL, $3);}
               ;
 
+forecast: FORECAST ';' {driver.forecast();}
+          | FORECAST '(' forecast_options ')' ';' {driver.forecast();}
+          | FORECAST symbol_list ';' {driver.forecast();}
+          | FORECAST '(' forecast_options ')' symbol_list ';' {driver.forecast();}
+          ;
+
+forecast_options: forecast_option
+          | forecast_options COMMA forecast_option
+          ;
+
+forecast_option: o_periods
+          | o_conf_sig
+          ;
+
 number : INT_NUMBER
        | FLOAT_NUMBER
        ;

preprocessor/DynareFlex.ll

@@ -130,6 +130,7 @@ int sigma_e = 0;
 <INITIAL>bvar_forecast {BEGIN DYNARE_STATEMENT; return token::BVAR_FORECAST; }
 <INITIAL>dynare_sensitivity {BEGIN DYNARE_STATEMENT; return token::DYNARE_SENSITIVITY;}
 <INITIAL>initval_file {BEGIN DYNARE_STATEMENT; return token::INITVAL_FILE;}
+<INITIAL>forecast {BEGIN DYNARE_STATEMENT; return token::FORECAST;}
  /* End of a Dynare statement */
 
 <DYNARE_STATEMENT>; {

preprocessor/ParsingDriver.cc

@@ -350,6 +350,14 @@ ParsingDriver::homotopy_val(string *name, NodeID val1, NodeID val2)
   delete name;
 }
 
+void
+ParsingDriver::forecast()
+{
+  mod_file->addStatement(new ForecastStatement(symbol_list, options_list));
+  symbol_list.clear();
+  options_list.clear();
+}
+
 void
 ParsingDriver::use_dll()
 {

preprocessor/include/ComputingTasks.hh

@@ -88,6 +88,18 @@ public:
   virtual void writeOutput(ostream &output, const string &basename) const;
 };
 
+class ForecastStatement : public Statement
+{
+private:
+  const SymbolList symbol_list;
+  const OptionsList options_list;
+public:
+  ForecastStatement(const SymbolList &symbol_list_arg,
+                      const OptionsList &options_list_arg);
+  virtual void checkPass(ModFileStructure &mod_file_struct);
+  virtual void writeOutput(ostream &output, const string &basename) const;
+};
+
 class RamseyPolicyStatement : public Statement
 {
 private:

preprocessor/include/ParsingDriver.hh

@@ -304,6 +304,8 @@ public:
   void optim_options_num(string *name, string *value);
   //! Prints varops instructions
   void set_varobs();
+  //! Forecast Statement
+  void forecast();
   void set_trends();
   void set_trend_element(string *arg1, NodeID arg2);
   void set_unit_root_vars();