1function [forcs, e]= mcforecast3(cL,H,mcValue,shocks,forcs,T,R,mv,mu)
2% [forcs, e] = mcforecast3(cL,H,mcValue,shocks,forcs,T,R,mv,mu)
3% Computes the shock values for constrained forecasts necessary to keep
4% endogenous variables at their constrained paths
5%
6% INPUTS
7%  o cL             [scalar]                            number of controlled periods
8%  o H              [scalar]                            number of forecast periods
9%  o mcValue        [n_controlled_vars by cL double]    paths for constrained variables
10%  o shocks         [nexo by H double]                  shock values draws (with zeros for controlled_varexo)
11%  o forcs          [n_endovars by H+1 double]          matrix of endogenous variables storing the inital condition
12%  o T              [n_endovars by n_endovars double]   transition matrix of the state equation.
13%  o R              [n_endovars by n_exo double]        matrix relating the endogenous variables to the innovations in the state equation.
14%  o mv             [n_controlled_exo by n_endovars boolean]   indicator vector  selecting constrained endogenous variables
15%  o mu             [n_controlled_vars by nexo boolean]        indicator vector selecting controlled exogenous variables
16% OUTPUTS
17%  o forcs          [n_endovars by H+1 double]          matrix of forecasted endogenous variables
18%  o e              [nexo by H double]                  matrix of exogenous variables
19%
20% Algorithm:
21%   Relies on state-space form:
22%       y_t=T*y_{t-1}+R*shocks(:,t)
23%   Shocks are split up into shocks_uncontrolled and shockscontrolled while
24%   the endogenous variables are also split up into controlled and
25%   uncontrolled ones to get:
26%       y_t(controlled_vars_index)=T*y_{t-1}(controlled_vars_index)+R(controlled_vars_index,uncontrolled_shocks_index)*shocks_uncontrolled_t
27%                    + R(controlled_vars_index,controlled_shocks_index)*shocks_controlled_t
28%
29%   This is then solved to get:
30%       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)
31%
32%   Variable number of controlled vars are allowed in different
33%   periods. Missing control information are indicated by NaN in
34%   y_t(controlled_vars_index).
35%
36%   After obtaining the shocks, and for uncontrolled periods, the state-space representation
37%       y_t=T*y_{t-1}+R*shocks(:,t)
38%   is used for forecasting
39%
40% Copyright (C) 2006-2017 Dynare Team
41%
42% This file is part of Dynare.
43%
44% Dynare is free software: you can redistribute it and/or modify
45% it under the terms of the GNU General Public License as published by
46% the Free Software Foundation, either version 3 of the License, or
47% (at your option) any later version.
48%
49% Dynare is distributed in the hope that it will be useful,
50% but WITHOUT ANY WARRANTY; without even the implied warranty of
51% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
52% GNU General Public License for more details.
53%
54% You should have received a copy of the GNU General Public License
55% along with Dynare.  If not, see <http://www.gnu.org/licenses/>.
56
57if cL
58    e = zeros(size(mcValue,1),cL);
59    for t = 1:cL
60        % missing conditional values are indicated by NaN
61        k = find(isfinite(mcValue(:,t)));
62        e(k,t) = inv(mv(k,:)*R*mu(:,k))*(mcValue(k,t)-mv(k,:)*T*forcs(:,t)-mv(k,:)*R*shocks(:,t));
63        forcs(:,t+1) = T*forcs(:,t)+R*(mu(:,k)*e(k,t)+shocks(:,t));
64    end
65end
66for t = cL+1:H
67    forcs(:,t+1) = T*forcs(:,t)+R*shocks(:,t);
68end