1function [yf,int_width,int_width_ME]=forcst(dr,y0,horizon,var_list,M_,oo_,options_)
2% function [yf,int_width,int_width_ME]=forecst(dr,y0,horizon,var_list,M_,oo_,options_)
3%   computes mean forecast for a given value of the parameters
4%   computes also confidence band for the forecast
5%
6% INPUTS:
7%   dr:          structure containing decision rules
8%   y0:          initial values
9%   horizon:     nbr of periods to forecast
10%   var_list:    list of variables (character matrix)
11%   M_:          Dynare model structure
12%   options_:    Dynare options structure
13%   oo_:         Dynare results structure
14
15% OUTPUTS:
16%   yf:          mean forecast
17%   int_width:   distance between upper bound and
18%                mean forecast
19%   int_width_ME:distance between upper bound and
20%                mean forecast when considering measurement error
21%
22% SPECIAL REQUIREMENTS
23%    none
24
25% Copyright (C) 2003-2019 Dynare Team
26%
27% This file is part of Dynare.
28%
29% Dynare is free software: you can redistribute it and/or modify
30% it under the terms of the GNU General Public License as published by
31% the Free Software Foundation, either version 3 of the License, or
32% (at your option) any later version.
33%
34% Dynare is distributed in the hope that it will be useful,
35% but WITHOUT ANY WARRANTY; without even the implied warranty of
36% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
37% GNU General Public License for more details.
38%
39% You should have received a copy of the GNU General Public License
40% along with Dynare.  If not, see <http://www.gnu.org/licenses/>.
41
42oo_=make_ex_(M_,options_,oo_);
43yf = simult_(M_,options_,y0,dr,zeros(horizon,M_.exo_nbr),1);
44nstatic = M_.nstatic;
45nspred = M_.nspred;
46nc = size(dr.ghx,2);
47endo_nbr = M_.endo_nbr;
48inv_order_var = dr.inv_order_var;
49[A,B] = kalman_transition_matrix(dr,nstatic+(1:nspred),1:nc,M_.exo_nbr);
50
51if isempty(var_list)
52    var_list = M_.endo_names(1:M_.orig_endo_nbr);
53end
54nvar = length(var_list);
55ivar = zeros(nvar, 1);
56for i=1:nvar
57    i_tmp = strmatch(var_list{i}, M_.endo_names, 'exact');
58    if isempty(i_tmp)
59        disp(var_list{i});
60        error ('One of the variable specified does not exist') ;
61    else
62        ivar(i) = i_tmp;
63    end
64end
65
66ghx1 = dr.ghx(inv_order_var(ivar),:);
67ghu1 = dr.ghu(inv_order_var(ivar),:);
68
69%initialize recursion
70sigma_u = B*M_.Sigma_e*B';
71sigma_u1 = ghu1*M_.Sigma_e*ghu1';
72sigma_y = 0; %no uncertainty about the states
73
74var_yf = NaN(horizon,nvar); %initialize
75for i = 1:horizon
76    %map uncertainty about states into uncertainty about observables
77    sigma_y1 = ghx1*sigma_y*ghx1'+sigma_u1;
78    var_yf(i,:) = diag(sigma_y1)';
79    if i == horizon
80        break
81    end
82    %update uncertainty about states
83    sigma_u = A*sigma_u*A';
84    sigma_y = sigma_y+sigma_u;
85end
86if nargout==3
87    var_yf_ME=var_yf;
88    [loc_H, loc_varlist] = ismember(options_.varobs', options_.varlist);
89    loc_varlist(loc_varlist==0) = [];
90    if ~isempty(loc_varlist)
91        var_yf_ME(:,loc_varlist) = var_yf(:,loc_varlist)+repmat(diag(M_.H(loc_H,loc_H))', horizon, 1);
92    end
93    int_width_ME = zeros(horizon, nvar);
94end
95
96fact = norminv((1-options_.forecasts.conf_sig)/2, 0, 1);
97
98int_width = zeros(horizon, nvar);
99for i = 1:nvar
100    int_width(:,i) = -fact*sqrt(var_yf(:,i));
101    if nargout==3
102        int_width_ME(:,i) = -fact*sqrt(var_yf_ME(:,i));
103    end
104end
105
106yf = yf(ivar,:);
107