diff --git a/doc/dynare.texi b/doc/dynare.texi
index 19d1c4e1de01c40e822b28f06f0fc4fa128a7822..2b558d3f93af2a1eb72e68f879fdcca92aab98b7 100644
--- a/doc/dynare.texi
+++ b/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
diff --git a/matlab/imcforecast.m b/matlab/imcforecast.m
index 17587ade11086e39c786b1f97094e8dd94d569a1..3c605b51c2039572917e6c52e823adeacc36b39c 100644
--- a/matlab/imcforecast.m
+++ b/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;
diff --git a/matlab/mcforecast3.m b/matlab/mcforecast3.m
index f3ce87c68015c6db32b9b78860dbc9788b186fd8..3e1de2888f89cc03e7e97460a396bc3323e38056 100644
--- a/matlab/mcforecast3.m
+++ b/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
diff --git a/tests/Makefile.am b/tests/Makefile.am
index 1eefe3098c30b9578abec5a08684c77fd7a66af1..e86dfb5e99511627387eca38a0c98079d01e5030 100644
--- a/tests/Makefile.am
+++ b/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 \
diff --git a/tests/conditional_forecasts/fs2000_conditional_forecast_initval.mod b/tests/conditional_forecasts/fs2000_conditional_forecast_initval.mod
new file mode 100644
index 0000000000000000000000000000000000000000..0d32149301ca257d6d6c8b1bdca13d04356d5952
--- /dev/null
+++ b/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