diff --git a/doc/manual/source/the-model-file.rst b/doc/manual/source/the-model-file.rst
index 8ee7c26e03a368c78903c3f94780db3f033c3839..4f7e6cf85cfa68b0bba0d6f649c8891a77338e2a 100644
--- a/doc/manual/source/the-model-file.rst
+++ b/doc/manual/source/the-model-file.rst
@@ -9907,7 +9907,7 @@ the :comm:`bvar_forecast` command.
variables. This is done using the reduced form first order
state-space representation of the DSGE model by finding the
structural shocks that are needed to match the restricted
- paths. Consider the an augmented state space representation that
+ paths. Consider the augmented state space representation that
stacks both predetermined and non-predetermined variables into a
vector :math:`y_{t}`:
@@ -9915,43 +9915,83 @@ the :comm:`bvar_forecast` command.
y_t=Ty_{t-1}+R\varepsilon_t
- Both :math:`y_t` and :math:`\varepsilon_t` are split up into
- controlled and uncontrolled ones to get:
+ Both :math:`y_t` and :math:`\varepsilon_t` are split up into controlled and
+ uncontrolled ones, and we assume without loss of generality that the
+ constrained endogenous variables and the controlled shocks come first :
.. math::
- y_t(contr\_vars)=Ty_{t-1}(contr\_vars)+R(contr\_vars,uncontr\_shocks)\varepsilon_t(uncontr\_shocks) \\
- + R(contr\_vars,contr\_shocks)\varepsilon_t(contr\_shocks)
-
- which can be solved algebraically for :math:`\varepsilon_t(contr\_shocks)`.
-
- Using these controlled shocks, the state-space representation can
- be used for forecasting. A few things need to be noted. First, it
- is assumed that controlled exogenous variables are fully under
- control of the policy maker for all forecast periods and not just
- for the periods where the endogenous variables are controlled. For
- all uncontrolled periods, the controlled exogenous variables are
- assumed to be 0. This implies that there is no forecast
- uncertainty arising from these exogenous variables in uncontrolled
- periods. Second, by making use of the first order state space
- solution, even if a higher-order approximation was performed, the
- conditional forecasts will be based on a first order
- approximation. Third, although controlled exogenous variables are
- taken as instruments perfectly under the control of the
- policy-maker, they are nevertheless random and unforeseen shocks
- from the perspective of the households. That is, households are in
- each period surprised by the realization of a shock that keeps the
- controlled endogenous variables at their respective level. Fourth,
- keep in mind that if the structural innovations are correlated,
- because the calibrated or estimated covariance matrix has non zero
- off diagonal elements, the results of the conditional forecasts
- will depend on the ordering of the innovations (as declared after
- ``varexo``). As in VAR models, a Cholesky decomposition is used to
- factorize the covariance matrix and identify orthogonal
- impulses. It is preferable to declare the correlations in the
- model block (explicitly imposing the identification restrictions),
- unless you are satisfied with the implicit identification
- restrictions implied by the Cholesky decomposition.
+ \begin{pmatrix}
+ y_{c,t}\\
+ y_{u,t}
+ \end{pmatrix}
+ =
+ \begin{pmatrix}
+ T_{c,c} & T_{c,u}\\
+ T_{u,c} & T_{u,u}
+ \end{pmatrix}
+ \begin{pmatrix}
+ y_{c,t-1}\\
+ y_{u,t-1}
+ \end{pmatrix}
+ +
+ \begin{pmatrix}
+ R_{c,c} & R_{c,u}\\
+ R_{u,c} & R_{u,u}
+ \end{pmatrix}
+ \begin{pmatrix}
+ \varepsilon_{c,t}\\
+ \varepsilon_{u,t}
+ \end{pmatrix}
+
+ where matrices :math:`T` and :math:`R` are partitioned consistently with the
+ vectors of endogenous variables and innovations. Provided that matrix
+ :math:`R_{c,c}` is square and full rank (a necessary condition is that the
+ number of free endogenous variables matches the number of free innovations),
+ given :math:`y_{c,t}`, :math:`\varepsilon_{u,t}` and :math:`y_{t-1}` the
+ first block of equations can be solved for :math:`\varepsilon_{c,t}`:
+
+ .. math::
+
+ \varepsilon_{c,t} = R_{c,c}^{-1}\bigl( y_{c,t} - T_{c,c}y_{c,t} - T_{c,u}y_{u,t} - R_{c,u}\varepsilon_{u,t}\bigr)
+
+ and :math:`y_{u,t}` can be updated by evaluating the second block of equations:
+
+ .. math::
+
+ y_{u,t} = T_{u,c}y_{c,t-1} + T_{u,u}y_{u,t-1} + R_{u,c}\varepsilon_{c,t} + R_{u,u}\varepsilon_{u,t}
+
+ By iterating over these two blocks of equations, we can build a forecast for
+ all the endogenous variables in the system conditional on paths for a subset of the
+ endogenous variables. If the distribution of the free innovations
+ :math:`\varepsilon_{u,t}` is provided (*i.e.* some of them have positive
+ variances) this exercise is replicated (the number of replication is
+ controlled by the option :opt:`replic` described below) by drawing different
+ sequences of free innovations. The result is a predictive distribution for
+ the uncontrolled endogenous variables, :math:`y_{u,t}`, that Dynare will use to report
+ confidence bands around the point conditional forecast.
+
+ A few things need to be noted. First, the controlled
+ exogenous variables are set to zero for the uncontrolled periods. This implies
+ that there is no forecast uncertainty arising from these exogenous variables
+ in uncontrolled periods. Second, by making use of the first order state
+ space solution, even if a higher-order approximation was performed, the
+ conditional forecasts will be based on a first order approximation. Since
+ the controlled exogenous variables are identified on the basis of the
+ reduced form model (*i.e.* after solving for the expectations), they are
+ unforeseen shocks from the perspective of the agents in the model. That is,
+ agents expect the endogenous variables to return to their respective steady
+ state levels but are surprised in each period by the realisation of shocks
+ keeping the endogenous variables along a predefined (unexpected) path.
+ Fourth, if the structural innovations are correlated, because the calibrated
+ or estimated covariance matrix has non zero off diagonal elements, the
+ results of the conditional forecasts will depend on the ordering of the
+ innovations (as declared after ``varexo``). As in VAR models, a Cholesky
+ decomposition is used to factorise the covariance matrix and identify
+ orthogonal impulses. It is preferable to declare the correlations in the
+ model block (explicitly imposing the identification restrictions), unless
+ you are satisfied with the implicit identification restrictions implied by
+ the Cholesky decomposition.
This command has to be called after ``estimation`` or ``stoch_simul``.
@@ -9982,7 +10022,7 @@ the :comm:`bvar_forecast` command.
.. option:: replic = INTEGER
- Number of simulations. Default: ``5000``.
+ Number of simulations used to compute the conditional forecast uncertainty. Default: ``5000``.
.. option:: conf_sig = DOUBLE
diff --git a/matlab/mcforecast3.m b/matlab/mcforecast3.m
index 47a0eab87c0cdc5aa6f3015190c17a250ce1217a..f249760e26b2b30ee5ed55b915593b2571d9c00d 100644
--- a/matlab/mcforecast3.m
+++ b/matlab/mcforecast3.m
@@ -1,43 +1,57 @@
-function [forcs, e]= mcforecast3(cL,H,mcValue,shocks,forcs,T,R,mv,mu)
-% [forcs, e] = mcforecast3(cL,H,mcValue,shocks,forcs,T,R,mv,mu)
+function [forcs, e] = 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
%
-% INPUTS
-% o cL [scalar] number of controlled periods
-% o H [scalar] number of forecast periods
-% o mcValue [n_controlled_vars by cL double] paths for constrained variables
-% o shocks [nexo by H double] shock values draws (with zeros for controlled_varexo)
-% o forcs [n_endovars by H+1 double] matrix of endogenous variables storing the inital condition
-% o T [n_endovars by n_endovars double] transition matrix of the state equation.
-% o R [n_endovars by n_exo double] matrix relating the endogenous variables to the innovations in the state equation.
-% o mv [n_controlled_exo by n_endovars boolean] indicator vector selecting constrained endogenous variables
-% o mu [n_controlled_vars by nexo boolean] indicator vector selecting controlled exogenous variables
-% OUTPUTS
-% o forcs [n_endovars by H+1 double] matrix of forecasted endogenous variables
-% o e [nexo by H double] matrix of exogenous variables
-%
-% Algorithm:
+% INPUTS:
+% - cL [integer] scalar, number of controlled periods
+% - H [integer] scalar, number of forecast periods
+% - mcValue [double] n_controlled_vars*cL array, paths for constrained variables
+% - shocks [double] n_controlled_vars*cL array, shock values draws (with zeros for controlled_varexo)
+% - forcs [double] n_endovars*(H+1) matrix of endogenous variables storing the inital condition
+% - T [double] n_endovars*n_endovars array, transition matrix of the state equation.
+% - R [double] n_endovars*n_exo array, matrix relating the endogenous variables to the innovations in the state equation.
+% - mv [logical] n_controlled_exo*n_endovars array, indicator selecting constrained endogenous variables
+% - mu [logical] n_controlled_vars*nexo array, indicator selecting controlled exogenous variables
+%
+% OUTPUTS:
+% - forcs [double] n_endovars*(H+1) array, forecasted endogenous variables
+% - e [double] nexo*H array, exogenous variables
+%
+% ALGORITHM:
+%
% Relies on state-space form:
-% y_t=T*y_{t-1}+R*shocks(:,t)
-% Shocks are split up into shocks_uncontrolled and shockscontrolled while
-% the endogenous variables are also split up into controlled and
-% uncontrolled ones to get:
-% y_t(controlled_vars_index)=T*y_{t-1}(controlled_vars_index)+R(controlled_vars_index,uncontrolled_shocks_index)*shocks_uncontrolled_t
-% + R(controlled_vars_index,controlled_shocks_index)*shocks_controlled_t
-%
-% This is then solved to get:
-% shocks_controlled_t=(y_t(controlled_vars_index)-(T*y_{t-1}(controlled_vars_index)+R(controlled_vars_index,uncontrolled_shocks_index)*shocks_uncontrolled_t)/R(controlled_vars_index,controlled_shocks_index)
-%
-% Variable number of controlled vars are allowed in different
-% periods. Missing control information are indicated by NaN in
-% y_t(controlled_vars_index).
-%
-% After obtaining the shocks, and for uncontrolled periods, the state-space representation
-% y_t=T*y_{t-1}+R*shocks(:,t)
-% is used for forecasting
%
-% Copyright (C) 2006-2017 Dynare Team
+% yₜ = T yₜ₋₁ + R εₜ
+%
+% Both yₜ, the vector of endogenous variables, and εₜ are split up into controlled
+% and uncontrolled ones, and we assume, without loss of generality, that the
+% constrained endogenous variables and the controlled shocks come first :
+%
+% ⎧ y₁ₜ ⎫ ⎧ T₁₁ T₁₂ ⎫ ⎧ y₁ₜ₋₁ ⎫ ⎧ R₁₁ R₁₂ ⎫ ⎧ ε₁ₜ ⎫
+% ⎩ y₂ₜ ⎭ = ⎩ T₂₁ T₂₂ ⎭ ⎩ y₂ₜ₋₁ ⎭ + ⎩ R₂₁ R₂₂ ⎭ ⎩ ε₂ₜ ⎭
+%
+% where matrices T and R are partitioned consistently with the
+% vectors of endogenous variables and innovations. Provided that matrix
+% R₁₁ is square and full rank (a necessary condition is that the
+% number of free endogenous variables matches the number of free innovations),
+% given y₁ₜ, ε₂ₜ and yₜ₋₁ the first block of equations can be solved for ε₁ₜ:
+%
+% ε₁ₜ = R₁₁ \ ( y₁ₜ - T₁₁y₁ₜ₋₁ - T₁₂y₂ₜ₋₁ - R₁₂ε₂ₜ )
+%
+% and y₂ₜ can be updated by evaluating the second block of equations:
+%
+% y₂ₜ = T₂₁y₁ₜ₋₁ + T₂₂y₂ₜ₋₁ + R₂₁ε₁ₜ + R₂₂ε₂ₜ
+%
+% By iterating over these two blocks of equations, we can build a forecast for
+% all the endogenous variables in the system conditional on paths for a subset of the
+% endogenous variables. This exercise is replicated by drawing different
+% sequences of free innovations. The result is a predictive distribution for
+% the uncontrolled endogenous variables, y₂ₜ, that Dynare will use to report
+% confidence bands around the point conditional forecast.
+% is used for forecasting
+
+% Copyright © 2006-2022 Dynare Team
%
% This file is part of Dynare.
%
@@ -54,15 +68,14 @@ function [forcs, e]= mcforecast3(cL,H,mcValue,shocks,forcs,T,R,mv,mu)
% You should have received a copy of the GNU General Public License
% along with Dynare. If not, see <https://www.gnu.org/licenses/>.
-if cL
- e = zeros(size(mcValue,1),cL);
- for t = 1:cL
- % missing conditional values are indicated by NaN
- k = find(isfinite(mcValue(:,t)));
- e(k,t) = inv(mv(k,:)*R*mu(:,k))*(mcValue(k,t)-mv(k,:)*T*forcs(:,t)-mv(k,:)*R*shocks(:,t));
- forcs(:,t+1) = T*forcs(:,t)+R*(mu(:,k)*e(k,t)+shocks(:,t));
+ if cL
+ e = zeros(size(mcValue,1),cL);
+ for t = 1:cL % Loop over the two blocks of equations
+ k = find(isfinite(mcValue(:,t))); % missing conditional values are indicated by NaN
+ e(k,t) = inv(mv(k,:)*R*mu(:,k))*(mcValue(k,t)-mv(k,:)*T*forcs(:,t)-mv(k,:)*R*shocks(:,t));
+ forcs(:,t+1) = T*forcs(:,t)+R*(mu(:,k)*e(k,t)+shocks(:,t));
+ end
end
-end
-for t = cL+1:H
- forcs(:,t+1) = T*forcs(:,t)+R*shocks(:,t);
-end
\ No newline at end of file
+ for t = cL+1:H
+ forcs(:,t+1) = T*forcs(:,t)+R*shocks(:,t);
+ end
\ No newline at end of file