Merge pull request #692 from JohannesPfeifer/imcforecast

Account for initial values when making conditional forecasts
time-shift
MichelJuillard 2015-07-20 10:18:43 +02:00
commit ce3dc068c3
8 changed files with 298 additions and 28 deletions

View File

@ -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

View File

@ -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

View File

@ -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;

View File

@ -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

View File

@ -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);

View File

@ -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 \

View File

@ -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

View File

@ -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