added possibility to do conditional forecast with a calibrated model + bug...

added possibility to do conditional forecast with a calibrated model + bug correction + new test file

doc/dynare.texi

@@ -4322,7 +4322,7 @@ Use @code{plot_conditional_forecast} to graph the results.
 
 @table @code
 
-@item parameter_set = @code{prior_mode} | @code{prior_mean} | @code{posterior_mode} | @code{posterior_mean} | @code{posterior_median}
+@item parameter_set = @code{calibration} | @code{prior_mode} | @code{prior_mean} | @code{posterior_mode} | @code{posterior_mean} | @code{posterior_median}
 Specify the parameter set to use for the forecasting. No default
 value, mandatory option.
 
@@ -4362,7 +4362,7 @@ end;
 
 conditional_forecast(parameter_set = calibration, controlled_varexo = (e, u), replic = 3000);
 
-plot_conditional_forecast(periods = 10) e u;
+plot_conditional_forecast(periods = 10) a y;
 @end example
 
 @end deffn

matlab/imcforecast.m

@@ -45,38 +45,10 @@ function imcforecast(constrained_paths, constrained_vars, options_cond_fcst)
 
 global options_ oo_ M_ bayestopt_
 
-if isfield(options_cond_fcst,'parameter_set') || isempty(options_cond_fcst.parameter_set)
+if ~isfield(options_cond_fcst,'parameter_set') || isempty(options_cond_fcst.parameter_set)
     options_cond_fcst.parameter_set = 'posterior_mode';
 end
 
-if ischar(options_cond_fcst.parameter_set)
-    switch options_cond_fcst.parameter_set
-      case 'posterior_mode'
-        xparam = get_posterior_parameters('mode');
-      case 'posterior_mean'
-        xparam = get_posterior_parameters('mean');
-      case 'posterior_median'
-        xparam = get_posterior_parameters('median');
-      case 'prior_mode'
-        xparam = bayestopt_.p5(:);
-      case 'prior_mean'
-        xparam = bayestopt_.p1;
-      otherwise
-        disp('eval_likelihood:: If the input argument is a string, then it has to be equal to:')
-        disp('                   ''posterior_mode'', ')
-        disp('                   ''posterior_mean'', ')
-        disp('                   ''posterior_median'', ')
-        disp('                   ''prior_mode'' or')
-        disp('                   ''prior_mean''.')
-        error('imcforecast:: Wrong argument type!')
-    end
-else
-    xparam = options_cond_fcst.parameter_set;
-    if length(xparam)~=length(M_.params)
-        error('imcforecast:: The dimension of the vector of parameters doesn''t match the number of estimated parameters!')
-    end
-end
-
 if ~isfield(options_cond_fcst,'replic') || isempty(options_cond_fcst.replic)
     options_cond_fcst.replic = 5000;
 end
@@ -89,56 +61,100 @@ if ~isfield(options_cond_fcst,'conf_sig') || isempty(options_cond_fcst.conf_sig)
     options_cond_fcst.conf_sig = .8;
 end
 
-set_parameters(xparam);
-
-n_varobs = size(options_.varobs,1);
-rawdata = read_variables(options_.datafile,options_.varobs,[],options_.xls_sheet,options_.xls_range);
-options_ = set_default_option(options_,'nobs',size(rawdata,1)-options_.first_obs+1);
-gend = options_.nobs;
-rawdata = rawdata(options_.first_obs:options_.first_obs+gend-1,:);
-% Transform the data.
-if options_.loglinear
-    if ~options_.logdata
-        rawdata = log(rawdata);  
-    end
-end
-% Test if the data set is real.
-if ~isreal(rawdata)
-    error('There are complex values in the data! Probably  a wrong transformation')
-end
-% Detrend the data.
-options_.missing_data = any(any(isnan(rawdata)));
-if options_.prefilter == 1
-    if options_.missing_data
-        bayestopt_.mean_varobs = zeros(n_varobs,1);
-        for variable=1:n_varobs
-            rdx = find(~isnan(rawdata(:,variable)));
-            m = mean(rawdata(rdx,variable));
-            rawdata(rdx,variable) = rawdata(rdx,variable)-m;
-            bayestopt_.mean_varobs(variable) = m;
+if isequal(options_cond_fcst.parameter_set,'calibration')
+    estimated_model = 0;
+else
+    estimated_model = 1;
+end
+
+if estimated_model
+    if ischar(options_cond_fcst.parameter_set)
+        switch options_cond_fcst.parameter_set
+          case 'posterior_mode'
+            xparam = get_posterior_parameters('mode');
+          case 'posterior_mean'
+            xparam = get_posterior_parameters('mean');
+          case 'posterior_median'
+            xparam = get_posterior_parameters('median');
+          case 'prior_mode'
+            xparam = bayestopt_.p5(:);
+          case 'prior_mean'
+            xparam = bayestopt_.p1;
+          otherwise
+            disp('imcforecast:: If the input argument is a string, then it has to be equal to:')
+            disp('                   ''calibration'', ')
+            disp('                   ''posterior_mode'', ')
+            disp('                   ''posterior_mean'', ')
+            disp('                   ''posterior_median'', ')
+            disp('                   ''prior_mode'' or')
+            disp('                   ''prior_mean''.')
+            error('imcforecast:: Wrong argument type!')
         end
     else
-        bayestopt_.mean_varobs = mean(rawdata,1)';
-        rawdata = rawdata-repmat(bayestopt_.mean_varobs',gend,1);
+        xparam = options_cond_fcst.parameter_set;
+        if length(xparam)~=length(M_.params)
+            error('imcforecast:: The dimension of the vector of parameters doesn''t match the number of estimated parameters!')
+        end
     end
-end
-data = transpose(rawdata);
-% Handle the missing observations.
-[data_index,number_of_observations,no_more_missing_observations] = describe_missing_data(data,gend,n_varobs);
-missing_value = ~(number_of_observations == gend*n_varobs);
-
-[atT,innov,measurement_error,filtered_state_vector,ys,trend_coeff] = DsgeSmoother(xparam,gend,data,data_index,number_of_observations);
 
-trend = repmat(ys,1,options_cond_fcst.periods+1);
-for i=1:M_.endo_nbr
-    j = strmatch(deblank(M_.endo_names(i,:)),options_.varobs,'exact');
-    if ~isempty(j)
-        trend(i,:) = trend(i,:)+trend_coeff(j)*(gend+(0:options_cond_fcst.periods));
+    set_parameters(xparam);
+
+    n_varobs = size(options_.varobs,1);
+    rawdata = read_variables(options_.datafile,options_.varobs,[],options_.xls_sheet,options_.xls_range);
+    options_ = set_default_option(options_,'nobs',size(rawdata,1)-options_.first_obs+1);
+    gend = options_.nobs;
+    rawdata = rawdata(options_.first_obs:options_.first_obs+gend-1,:);
+    % Transform the data.
+    if options_.loglinear
+        if ~options_.logdata
+            rawdata = log(rawdata);  
+        end
+    end
+    % Test if the data set is real.
+    if ~isreal(rawdata)
+        error('There are complex values in the data! Probably  a wrong transformation')
+    end
+    % Detrend the data.
+    options_.missing_data = any(any(isnan(rawdata)));
+    if options_.prefilter == 1
+        if options_.missing_data
+            bayestopt_.mean_varobs = zeros(n_varobs,1);
+            for variable=1:n_varobs
+                rdx = find(~isnan(rawdata(:,variable)));
+                m = mean(rawdata(rdx,variable));
+                rawdata(rdx,variable) = rawdata(rdx,variable)-m;
+                bayestopt_.mean_varobs(variable) = m;
+            end
+        else
+            bayestopt_.mean_varobs = mean(rawdata,1)';
+            rawdata = rawdata-repmat(bayestopt_.mean_varobs',gend,1);
+        end
+    end
+    data = transpose(rawdata);
+    % Handle the missing observations.
+    [data_index,number_of_observations,no_more_missing_observations] = describe_missing_data(data,gend,n_varobs);
+    missing_value = ~(number_of_observations == gend*n_varobs);
+
+    [atT,innov,measurement_error,filtered_state_vector,ys,trend_coeff] = DsgeSmoother(xparam,gend,data,data_index,number_of_observations);
+
+    trend = repmat(ys,1,options_cond_fcst.periods+1);
+    for i=1:M_.endo_nbr
+        j = strmatch(deblank(M_.endo_names(i,:)),options_.varobs,'exact');
+        if ~isempty(j)
+            trend(i,:) = trend(i,:)+trend_coeff(j)*(gend+(0:options_cond_fcst.periods));
+        end
     end
+    trend = trend(oo_.dr.order_var,:);
+
+    InitState(:,1) = atT(:,end);
+else
+    InitState(:,1) = zeros(M_.endo_nbr,1);
+    trend = repmat(oo_.steady_state(oo_.dr.order_var),1,options_cond_fcst.periods+1);
 end
-trend = trend(oo_.dr.order_var,:);
 
-InitState(:,1) = atT(:,end);
+if isempty(options_.qz_criterium)
+    options_.qz_criterium = 1+1e-6;
+end
 [T,R,ys,info] = dynare_resolve;
 
 sQ = sqrt(M_.Sigma_e);
@@ -146,9 +162,7 @@ sQ = sqrt(M_.Sigma_e);
 NumberOfStates = length(InitState);
 FORCS1 = zeros(NumberOfStates,options_cond_fcst.periods+1,options_cond_fcst.replic);
 
-for b=1:options_cond_fcst.replic
-    FORCS1(:,1,b) = InitState;
-end
+FORCS1(:,1,:) = repmat(InitState,1,options_cond_fcst.replic);
 
 EndoSize = M_.endo_nbr;
 ExoSize = M_.exo_nbr;
@@ -164,7 +178,7 @@ idx = [];
 jdx = [];
 
 for i = 1:n1
-    idx = [idx ; oo_.dr.inv_order_var(strmatch(deblank(constrained_vars(i,:)),M_.endo_names,'exact'))];
+    idx = [idx ; oo_.dr.inv_order_var(constrained_vars(i,:))];
     jdx = [jdx ; strmatch(deblank(options_cond_fcst.controlled_varexo(i,:)),M_.exo_names,'exact')];
 end
 mv = zeros(n1,NumberOfStates);
@@ -182,7 +196,7 @@ end
 
 constrained_paths = bsxfun(@minus,constrained_paths,trend(idx,1:cL));
 
-randn('state',0);
+%randn('state',0);
 
 for b=1:options_cond_fcst.replic
     shocks = sQ*randn(ExoSize,options_cond_fcst.periods);

preprocessor/DynareBison.yy

@@ -94,7 +94,7 @@ class ParsingDriver;
 %token BVAR_PRIOR_MU BVAR_PRIOR_OMEGA BVAR_PRIOR_TAU BVAR_PRIOR_TRAIN
 %token BVAR_REPLIC BYTECODE
 %token CHANGE_TYPE CHECK CONDITIONAL_FORECAST CONDITIONAL_FORECAST_PATHS CONF_SIG CONSTANT CONTROLLED_VAREXO CORR COVAR CUTOFF
-%token DATAFILE DR_ALGO DROP DSAMPLE DYNASAVE DYNATYPE
+%token DATAFILE DR_ALGO DROP DSAMPLE DYNASAVE DYNATYPE CALIBRATION
 %token END ENDVAL EQUAL ESTIMATION ESTIMATED_PARAMS ESTIMATED_PARAMS_BOUNDS ESTIMATED_PARAMS_INIT
 %token FILENAME FILTER_STEP_AHEAD FILTERED_VARS FIRST_OBS
 %token <string_val> FLOAT_NUMBER
@@ -1969,6 +1969,8 @@ o_parameter_set : PARAMETER_SET EQUAL PRIOR_MODE
                   { driver.option_str("parameter_set", "posterior_mode"); }
                 | PARAMETER_SET EQUAL POSTERIOR_MEDIAN
                   { driver.option_str("parameter_set", "posterior_median"); }
+                | PARAMETER_SET EQUAL CALIBRATION
+                  { driver.option_str("parameter_set", "calibration"); }
                 ;
 o_shocks : SHOCKS EQUAL '(' list_of_symbol_lists ')' { driver.option_symbol_list("shocks"); };
 o_labels : LABELS EQUAL '(' symbol_list ')' { driver.option_symbol_list("labels"); };

preprocessor/DynareFlex.ll

@@ -443,6 +443,7 @@ string eofbuff;
 <DYNARE_STATEMENT>mh_recover {return token::MH_RECOVER;}
 <DYNARE_STATEMENT>planner_discount {return token::PLANNER_DISCOUNT;}
 <DYNARE_STATEMENT>labels {return token::LABELS;}
+<DYNARE_STATEMENT>calibration {return token::CALIBRATION;}
 
 <DYNARE_BLOCK>equation {return token::EQUATION;}
 <DYNARE_BLOCK>exclusion {return token::EXCLUSION;}

preprocessor/Shocks.cc

@@ -289,9 +289,9 @@ ConditionalForecastPathsStatement::writeOutput(ostream &output, const string &ba
        it != paths.end(); it++)
     {
       if (it == paths.begin())
-        output << "constrained_vars_ = '" << it->first << "';" << endl;
+        output << "constrained_vars_ = " << it->first << ";" << endl;
       else
-        output << "constrained_vars_ = char(constrained_vars_, '" << it->first << "');" << endl;
+        output << "constrained_vars_ = [constrained_vars_; " << it->first << "];" << endl;
 
       const vector<AbstractShocksStatement::DetShockElement> &elems = it->second;
       for (int i = 0; i < (int) elems.size(); i++)

tests/conditional_forecasts/fs2000.mod

+// See fs2000.mod in the examples/ directory for details on the model
+
+var m P c e W R k d n l gy_obs gp_obs y dA;
+varexo e_a e_m;
+
+parameters alp bet gam mst rho psi del;
+
+alp = 0.33;
+bet = 0.99;
+gam = 0.003;
+mst = 1.011;
+rho = 0.7;
+psi = 0.787;
+del = 0.02;
+
+model;
+dA = exp(gam+e_a);
+log(m) = (1-rho)*log(mst) + rho*log(m(-1))+e_m;
+-P/(c(+1)*P(+1)*m)+bet*P(+1)*(alp*exp(-alp*(gam+log(e(+1))))*k^(alp-1)*n(+1)^(1-alp)+(1-del)*exp(-(gam+log(e(+1)))))/(c(+2)*P(+2)*m(+1))=0;
+W = l/n;
+-(psi/(1-psi))*(c*P/(1-n))+l/n = 0;
+R = P*(1-alp)*exp(-alp*(gam+e_a))*k(-1)^alp*n^(-alp)/W;
+1/(c*P)-bet*P*(1-alp)*exp(-alp*(gam+e_a))*k(-1)^alp*n^(1-alp)/(m*l*c(+1)*P(+1)) = 0;
+c+k = exp(-alp*(gam+e_a))*k(-1)^alp*n^(1-alp)+(1-del)*exp(-(gam+e_a))*k(-1);
+P*c = m;
+m-1+d = l;
+e = exp(e_a);
+y = k(-1)^alp*n^(1-alp)*exp(-alp*(gam+e_a));
+gy_obs = dA*y/y(-1);
+gp_obs = (P/P(-1))*m(-1)/dA;
+end;
+
+initval;
+k = 6;
+m = mst;
+P = 2.25;
+c = 0.45;
+e = 1;
+W = 4;
+R = 1.02;
+d = 0.85;
+n = 0.19;
+l = 0.86;
+y = 0.6;
+gy_obs = exp(gam);
+gp_obs = exp(-gam);
+dA = exp(gam);
+end;
+
+shocks;
+var e_a; stderr 0.014;
+var e_m; stderr 0.005;
+end;
+
+steady;
+
+check;
+
+stoch_simul(irf=0);
+
+conditional_forecast_paths;
+var gy_obs;
+periods  1  2  3:5;
+values   0.01 -0.02 0;
+var gp_obs;
+periods 1:5;
+values  0.05;
+end;
+
+conditional_forecast(parameter_set=calibration, controlled_varexo=(e_a,e_m));
+
+plot_conditional_forecast(periods=10) gy_obs gp_obs;
\ No newline at end of file