diff --git a/doc/dynare.texi b/doc/dynare.texi index 8d4d7f874..34c3948a9 100644 --- a/doc/dynare.texi +++ b/doc/dynare.texi @@ -2169,20 +2169,24 @@ commands (@code{stoch_simul}, @code{estimation}@dots{}). It is not necessary to declare @code{0} as initial value for exogenous stochastic variables, since it is the only possible value. -This steady state will be used as the initial condition at all the -periods preceeding the first simulation period for the two possible +The subsequently computed steady state (not the initial values, use @xref{histval} or this) will be used as the initial condition at all the +periods preceeding the first simulation period for the three possible types of simulations in stochastic mode: @itemize @item -in @code{stoch_simul}, if the @code{periods} options is specified +in @xref{stoch_simul}, if the @code{periods} options is specified @item -in @code{forecast} (in this case, note that it is still possible to -modify some of these initial values with @code{histval}) +in @xref{forecast} as the initial point at which the forecasts are computed + +@item +in @xref{conditional_forecast} as the initial point at which the conditional forecasts are computed @end itemize +To start simulations at a particular set of starting values that are not a computed steady state, use @xref{histval}. + @optionshead @table @code @@ -2395,16 +2399,30 @@ before (except when a @code{steady} command doesn't follow an @customhead{In a stochastic simulation context} In the context of stochastic simulations, @code{histval} allows setting -the starting point of those simulations in the state space (it does not -affect the starting point for impulse response functions). As for the case of +the starting point of those simulations in the state space. As for the case of perfect foresight simulations, all not explicitly specified variables are set to 0. -Moreover, as only states enter the recursive policy functions, all values specified for control variables will be ignored. +Moreover, as only states enter the recursive policy functions, all values specified for control variables will be ignored. This can be used -For @ref{Ramsey} policy, it also specifies the values of the endogenous states at +@itemize + +@item +in @ref{stoch_simul}, if the @code{periods} options is specified. Note that this only affects the starting point for the simulation, but not for the impulse response functions. + +@item +in @ref{forecast} as the initial point at which the forecasts are computed + +@item +in @ref{conditional_forecast} as the initial point at which the conditional forecasts are computed + +@item +in @ref{Ramsey} policy, where it also specifies the values of the endogenous states at which the objective function of the planner is computed. Note that the initial values of the Lagrange multipliers associated with the planner's problem cannot be set (@pxref{planner_objective_value}). +@end itemize + + @examplehead @example @@ -3600,7 +3618,7 @@ solution is computed using the @code{extended_path} command. @deffn Command stoch_simul [@var{VARIABLE_NAME}@dots{}]; @deffnx Command stoch_simul (@var{OPTIONS}@dots{}) [@var{VARIABLE_NAME}@dots{}]; - +@anchor{stoch_simul} @descriptionhead @code{stoch_simul} solves a stochastic (@i{i.e.} rational @@ -6317,7 +6335,7 @@ Fields are of the form: @deffn Command conditional_forecast (@var{OPTIONS}@dots{}) [@var{VARIABLE_NAME}@dots{}]; - +@anchor{conditional_forecast} @descriptionhead This command computes forecasts on an estimated or calibrated model for a @@ -6434,6 +6452,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 @@ -6473,8 +6500,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/dyn_forecast.m b/matlab/dyn_forecast.m index b6ad8c960..30b867476 100644 --- a/matlab/dyn_forecast.m +++ b/matlab/dyn_forecast.m @@ -64,7 +64,7 @@ switch task horizon = 5; end if isempty(M_.endo_histval) - y0 = repmat(oo_.steady_state,1,maximum_lag); + y0 = repmat(oo_.dr.ys,1,maximum_lag); else y0 = M_.endo_histval; end diff --git a/matlab/imcforecast.m b/matlab/imcforecast.m index 17587ade1..fece7c82f 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,26 @@ 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 + if isempty(M_.endo_histval) + y0 = ys; + else + y0 = M_.endo_histval; + end + InitState(:,1) = y0(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 +180,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 +208,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 f3ce87c68..3e1de2888 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/matlab/plot_icforecast.m b/matlab/plot_icforecast.m index 84db7297b..7e4b171c3 100644 --- a/matlab/plot_icforecast.m +++ b/matlab/plot_icforecast.m @@ -10,7 +10,7 @@ function plot_icforecast(Variables,periods,options_) % SPECIAL REQUIREMENTS % This routine has to be called after imcforecast.m. -% Copyright (C) 2006-2013 Dynare Team +% Copyright (C) 2006-2014 Dynare Team % % This file is part of Dynare. % @@ -34,8 +34,16 @@ end load conditional_forecasts; +eval(['forecast_periods = length(forecasts.cond.Mean.' Variables(1,:) ');']); if nargin==1 || isempty(periods) % Set default number of periods. - eval(['periods = length(forecasts.cond.Mean.' Variables(1,:) ');']); + periods=forecast_periods; +else + periods=periods+1; %take care of initial period that is not forecasted +end + +if periods>forecast_periods + fprintf('\nplot_icforecast:: Number of periods for plotting exceeds forecast horizon. Setting periods to %d.\n',forecast_periods-1) + periods=forecast_periods; end forecasts.graph.OutputDirectoryName = CheckPath('graphs',forecasts.graph.fname); diff --git a/tests/Makefile.am b/tests/Makefile.am index 96b17a5b8..569c67161 100644 --- a/tests/Makefile.am +++ b/tests/Makefile.am @@ -135,6 +135,8 @@ 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 \ + conditional_forecasts/fs2000_conditional_forecast_histval.mod \ recursive/ls2003.mod \ recursive/ls2003_bayesian.mod \ recursive/ls2003_bayesian_xls.mod \ diff --git a/tests/conditional_forecasts/fs2000_conditional_forecast_histval.mod b/tests/conditional_forecasts/fs2000_conditional_forecast_histval.mod new file mode 100644 index 000000000..a5c6c9dcd --- /dev/null +++ b/tests/conditional_forecasts/fs2000_conditional_forecast_histval.mod @@ -0,0 +1,105 @@ +/* Mod file tests the correctness of the conditional_forecast command when used together with histval 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 and all other states to steady state +histval; +k(0) = 6; +m(0)=1.01100000000000; +P(0)=2.25815456051727; +y(0)=0.580764879486831; +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-8 + error('Unconditional Forecasts do not match') +end + +%compare conditional forecasts; histval here sets initval condition for capital different from steady state +initial_condition_states = oo_.dr.ys; +initial_condition_states(strmatch('k',M_.endo_names,'exact')) = 6; +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-8 + error('Unconditional Forecasts do not match') +end 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 000000000..f37287e2a --- /dev/null +++ b/tests/conditional_forecasts/fs2000_conditional_forecast_initval.mod @@ -0,0 +1,104 @@ +/* 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 + * - initval should not play a role as initial conditions need to be set with histval; with initval + * the subsequently computed steady state is used + */ + +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-8 + error('Unconditional Forecasts do not match') +end + +%compare conditional forecasts; initval should not play a role as initial +%conditions need to be set with histval; +initial_condition_states = oo_.dr.ys; +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-8 + error('Unconditional Forecasts do not match') +end