Account for initial values when making conditional forecasts

- Closes #655
- Implements unit test
- Saves distribution of controlled exogenous variables
- Adds information that the full path of constrained variables must be specified

doc/dynare.texi

@@ -6364,6 +6364,15 @@ Variable set by the @code{conditional_forecast} command. Stores the names of the
 Variable set by the @code{conditional_forecast} command. Stores the position of the constrained endogenous variables in declaration order.
 @end defvr
 
+@defvr {MATLAB/Octave variable} forecasts.controlled_exo_variables
+Variable set by the @code{conditional_forecast} command. Stores the values of the controlled exogenous 
+variables underlying the conditional forecasts to achieve the constrained endogenous 
+variables. Fields are number of constrained periods by 1 vectors and are of the form:
+@example
+@code{forecasts.controlled_exo_variables.@var{FORECAST_MOMENT}.@var{SHOCK_NAME}}
+@end example
+@end defvr
+
 @defvr {MATLAB/Octave variable} forecasts.graphs
 Variable set by the @code{conditional_forecast} command. Stores the information for generating the conditional forecast plots.
 @end defvr
@@ -6403,8 +6412,12 @@ Describes the path of constrained endogenous, before calling
 shocks in @code{shocks}, see @code{conditional_forecast} for an
 example.
 
-The syntax of the block is the same than the deterministic shocks in
-the @code{shocks} blocks (@pxref{Shocks on exogenous variables}).
+The syntax of the block is the same as for the deterministic shocks in
+the @code{shocks} blocks (@pxref{Shocks on exogenous variables}). Note that you need to specify the full path for all constrained endogenous 
+variables between the first and last specified period. If an intermediate period 
+is not specified, a value of 0 is assumed. That is, if you specify only 
+values for periods 1 and 3, the values for period 2 will be 0. Currently, it is not 
+possible to have uncontrolled intermediate periods.
 
 
 @end deffn

matlab/imcforecast.m

@@ -118,8 +118,6 @@ if estimated_model
     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);
     graph_title='Calibration';
 end
 
@@ -132,13 +130,21 @@ end
 if ~isdiagonal(M_.Sigma_e)
     warning(sprintf('The innovations are correlated (the covariance matrix has non zero off diagonal elements), the results of the conditional forecasts will\ndepend on the ordering of the innovations (as declared after varexo) because a Cholesky decomposition is used to factorize the covariance matrix.\n\n=> It is preferable to declare the correlations in the model block (explicitly imposing the identification restrictions), unless you are satisfied\nwith the implicit identification restrictions implied by the Cholesky decomposition.'))
 end
-
 sQ = chol(M_.Sigma_e,'lower');
 
+if ~estimated_model
+    InitState(:,1) = oo_.steady_state(oo_.dr.order_var)-ys(oo_.dr.order_var,:); %initial state in deviations from steady state
+    trend = repmat(ys(oo_.dr.order_var,:),1,options_cond_fcst.periods+1); %trend needs to contain correct steady state
+end
+sQ = sqrt(M_.Sigma_e);
+
+
+
+
 NumberOfStates = length(InitState);
 FORCS1 = zeros(NumberOfStates,options_cond_fcst.periods+1,options_cond_fcst.replic);
 
-FORCS1(:,1,:) = repmat(InitState,1,options_cond_fcst.replic);
+FORCS1(:,1,:) = repmat(InitState,1,options_cond_fcst.replic); %set initial steady state to deviations from steady state in first period
 
 EndoSize = M_.endo_nbr;
 ExoSize = M_.exo_nbr;
@@ -169,15 +175,19 @@ cL = size(constrained_paths,2);
 
 constrained_paths = bsxfun(@minus,constrained_paths,trend(idx,1:cL));
 
+FORCS1_shocks = zeros(n1,cL,options_cond_fcst.replic);
+
 %randn('state',0);
 
-for b=1:options_cond_fcst.replic
+for b=1:options_cond_fcst.replic %conditional forecast using cL set to constrained values
     shocks = sQ*randn(ExoSize,options_cond_fcst.periods);
     shocks(jdx,:) = zeros(length(jdx),options_cond_fcst.periods);
-    FORCS1(:,:,b) = mcforecast3(cL,options_cond_fcst.periods,constrained_paths,shocks,FORCS1(:,:,b),T,R,mv, mu)+trend;
+    [FORCS1(:,:,b), FORCS1_shocks(:,:,b)] = mcforecast3(cL,options_cond_fcst.periods,constrained_paths,shocks,FORCS1(:,:,b),T,R,mv, mu);
+    FORCS1(:,:,b)=FORCS1(:,:,b)+trend; %add trend
 end
 
 mFORCS1 = mean(FORCS1,3);
+mFORCS1_shocks = mean(FORCS1_shocks,3);
 
 tt = (1-options_cond_fcst.conf_sig)/2;
 t1 = round(options_cond_fcst.replic*tt);
@@ -193,16 +203,21 @@ for i = 1:EndoSize
           ' = [tmp(t1,:)'' ,tmp(t2,:)'' ]'';']);
 end
 
-clear FORCS1;
+for i = 1:n1
+    eval(['forecasts.controlled_exo_variables.Mean.' deblank(options_cond_fcst.controlled_varexo(i,:)) ' = mFORCS1_shocks(i,:)'';']);
+    tmp = sort(squeeze(FORCS1_shocks(i,:,:))');
+    eval(['forecasts.controlled_exo_variables.ci.' deblank(options_cond_fcst.controlled_varexo(i,:))  ...
+          ' = [tmp(t1,:)'' ,tmp(t2,:)'' ]'';']);
+end
+
+clear FORCS1 mFORCS1_shocks;
 
 FORCS2 = zeros(NumberOfStates,options_cond_fcst.periods+1,options_cond_fcst.replic);
-for b=1:options_cond_fcst.replic
-    FORCS2(:,1,b) = InitState;
-end
+FORCS2(:,1,:) = repmat(InitState,1,options_cond_fcst.replic); %set initial steady state to deviations from steady state in first period
 
 %randn('state',0);
 
-for b=1:options_cond_fcst.replic
+for b=1:options_cond_fcst.replic %conditional forecast using cL set to 0
     shocks = sQ*randn(ExoSize,options_cond_fcst.periods);
     shocks(jdx,:) = zeros(length(jdx),options_cond_fcst.periods);
     FORCS2(:,:,b) = mcforecast3(0,options_cond_fcst.periods,constrained_paths,shocks,FORCS2(:,:,b),T,R,mv, mu)+trend;

matlab/mcforecast3.m

@@ -1,4 +1,4 @@
-function forcs = mcforecast3(cL,H,mcValue,shocks,forcs,T,R,mv,mu)
+function [forcs, e]= mcforecast3(cL,H,mcValue,shocks,forcs,T,R,mv,mu)
 % forcs = mcforecast3(cL,H,mcValue,shocks,forcs,T,R,mv,mu)
 % Computes the shock values for constrained forecasts necessary to keep
 % endogenous variables at their constrained paths

tests/Makefile.am

@@ -134,6 +134,7 @@ MODFILES = \
 	simul/simul_ZLB_purely_forward_no_solution.mod \
 	conditional_forecasts/fs2000_cal.mod \
 	conditional_forecasts/fs2000_est.mod \
+	conditional_forecasts/fs2000_conditional_forecast_initval.mod \
 	recursive/ls2003.mod \
 	recursive/ls2003_bayesian.mod \
 	recursive/ls2003_bayesian_xls.mod \

tests/conditional_forecasts/fs2000_conditional_forecast_initval.mod

@@ -0,0 +1,101 @@
+/* Mod file tests the correctness of the conditional_forecast command when used together with initval by
+ * - checking whether the unconditional forecast from the conditional_forecast-command 
+ *      coincides with the one from the forecast command
+ * - checking whether the conditional forecast coincides with the path of 
+ *      capital derived when simulating the model with simult_ and the computed exogenous instruments
+*/
+
+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  2  3:5;
+values   1 1.04 0.98;
+end;
+
+%set capital to non-steady state value
+initval;
+k = 6; 
+end;
+
+conditional_forecast(periods=100,parameter_set=calibration,replic=10000, controlled_varexo=(e_m,e_a));
+
+plot_conditional_forecast(periods=100) gy_obs gp_obs k;
+forecast(periods=100);
+
+%compare unconditional forecasts
+cond_forecast=load('conditional_forecasts.mat');
+if max(abs(cond_forecast.forecasts.uncond.Mean.k(2:end)-oo_.forecast.Mean.k))>1e-3
+    error('Unconditional Forecasts do not match')
+end
+        
+%compare conditional forecasts
+initial_condition_states = oo_.steady_state;
+shock_matrix = zeros(options_cond_fcst_.periods ,M_.exo_nbr); %create shock matrix with found controlled shocks
+shock_matrix(1:5,strmatch('e_a',M_.exo_names,'exact')) = cond_forecast.forecasts.controlled_exo_variables.Mean.e_a; %set controlled shocks to their values
+shock_matrix(1:5,strmatch('e_m',M_.exo_names,'exact')) = cond_forecast.forecasts.controlled_exo_variables.Mean.e_m; %set controlled shocks to their values
+
+y_simult = simult_(initial_condition_states,oo_.dr,shock_matrix,1);
+
+if max(abs(y_simult(strmatch('k',M_.endo_names,'exact'),:)'-cond_forecast.forecasts.cond.Mean.k))>1e-3
+    error('Unconditional Forecasts do not match')
+end