1function [y_,int_width,int_width_ME]=simultxdet(y0,ex,ex_det, iorder,var_list,M_,oo_,options_)
2%function [y_,int_width]=simultxdet(y0,ex,ex_det, iorder,var_list,M_,oo_,options_)
3%
4% Simulates a stochastic model in the presence of deterministic exogenous shocks
5%
6% INPUTS:
7%    y0:        initial values, of length M_.maximum_lag
8%    ex:        matrix of stochastic exogenous shocks, starting at period 1
9%    ex_det:    matrix of deterministic exogenous shocks, starting at period 1-M_.maximum_lag
10%    iorder:    order of approximation
11%    var_list:  list of endogenous variables to simulate
12%   int_width_ME:distance between upper bound and
13%                mean forecast when considering measurement error
14% OUTPUTS:
15%   yf:          mean forecast
16%   int_width:   distance between upper bound and
17%                mean forecast
18%   int_width_ME:distance between upper bound and
19%                mean forecast when considering measurement error
20%
21% The forecast horizon is equal to size(ex, 1).
22% The condition size(ex,1)+M_.maximum_lag=size(ex_det,1) must be verified
23%  for consistency.
24
25% Copyright (C) 2008-2018 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
42dr = oo_.dr;
43ykmin = M_.maximum_lag;
44endo_nbr = M_.endo_nbr;
45exo_det_steady_state = oo_.exo_det_steady_state;
46nstatic = M_.nstatic;
47nspred = M_.nspred;
48nc = size(dr.ghx,2);
49iter = size(ex,1);
50if size(ex_det, 1) ~= iter+ykmin
51    error('Size mismatch: number of forecasting periods for stochastic exogenous and deterministic exogenous don''t match')
52end
53nx = size(dr.ghu,2);
54y_ = zeros(size(y0,1),iter+ykmin);
55y_(:,1:ykmin) = y0;
56k1 = [ykmin:-1:1];
57k2 = dr.kstate(find(dr.kstate(:,2) <= ykmin+1),[1 2]);
58k2 = k2(:,1)+(ykmin+1-k2(:,2))*endo_nbr;
59k3 = M_.lead_lag_incidence(1:ykmin,:)';
60k3 = find(k3(:));
61k4 = dr.kstate(find(dr.kstate(:,2) < ykmin+1),[1 2]);
62k4 = k4(:,1)+(ykmin+1-k4(:,2))*endo_nbr;
63
64nvar = length(var_list);
65if nvar == 0
66    nvar = endo_nbr;
67    ivar = [1:nvar];
68else
69    ivar=zeros(nvar,1);
70    for i=1:nvar
71        i_tmp = strmatch(var_list{i}, M_.endo_names, 'exact');
72        if isempty(i_tmp)
73            disp(var_list{i})
74            error (['One of the variable specified does not exist']) ;
75        else
76            ivar(i) = i_tmp;
77        end
78    end
79end
80
81if iorder == 1
82    for i = ykmin+1: iter+ykmin
83        tempx1 = y_(dr.order_var,k1);
84        tempx2 = tempx1-repmat(dr.ys(dr.order_var),1,ykmin);
85        tempx = tempx2(k2);
86        y_(dr.order_var,i) = dr.ys(dr.order_var)+dr.ghx*tempx+dr.ghu*ex(i-ykmin,:)';
87        for j=1:min(ykmin+M_.exo_det_length+1-i,M_.exo_det_length)
88            y_(dr.order_var,i) = y_(dr.order_var,i) + dr.ghud{j}*(ex_det(i+j-1,:)'-exo_det_steady_state);
89        end
90
91        k1 = k1+1;
92    end
93elseif iorder == 2
94    for i = ykmin+1: iter+ykmin
95        tempx1 = y_(dr.order_var,k1);
96        tempx2 = tempx1-repmat(dr.ys(dr.order_var),1,ykmin);
97        tempx = tempx2(k2);
98        tempu = ex(i-ykmin,:)';
99        tempuu = kron(tempu,tempu);
100        tempxx = kron(tempx,tempx);
101        tempxu = kron(tempx,tempu);
102        y_(dr.order_var,i) = dr.ys(dr.order_var)+dr.ghs2/2+dr.ghx*tempx+ ...
103            dr.ghu*tempu+0.5*(dr.ghxx*tempxx+dr.ghuu*tempuu)+dr.ghxu* ...
104            tempxu;
105        for j=1:min(ykmin+M_.exo_det_length+1-i,M_.exo_det_length)
106            tempud = ex_det(i+j-1,:)'-exo_det_steady_state;
107            tempudud = kron(tempud,tempud);
108            tempxud = kron(tempx,tempud);
109            tempuud = kron(tempu,tempud);
110            y_(dr.order_var,i) = y_(dr.order_var,i) + dr.ghud{j}*tempud + ...
111                dr.ghxud{j}*tempxud + dr.ghuud{j}*tempuud + ...
112                0.5*dr.ghudud{j,j}*tempudud;
113            for k=1:j-1
114                tempudk = ex_det(i+k-1,:)'-exo_det_steady_state;
115                tempududk = kron(tempudk,tempud);
116                y_(dr.order_var,i) = y_(dr.order_var,i) + ...
117                    dr.ghudud{k,j}*tempududk;
118            end
119        end
120        k1 = k1+1;
121    end
122end
123
124[A,B] = kalman_transition_matrix(dr,nstatic+(1:nspred),1:nc,M_.exo_nbr);
125
126inv_order_var = dr.inv_order_var;
127ghx1 = dr.ghx(inv_order_var(ivar),:);
128ghu1 = dr.ghu(inv_order_var(ivar),:);
129
130sigma_u = B*M_.Sigma_e*B';
131sigma_u1 = ghu1*M_.Sigma_e*ghu1';
132sigma_y = 0;
133
134var_yf=NaN(iter,nvar); %initialize
135for i=1:iter
136    sigma_y1 = ghx1*sigma_y*ghx1'+sigma_u1;
137    var_yf(i,:) = diag(sigma_y1)';
138    if i == iter
139        break
140    end
141    sigma_u = A*sigma_u*A';
142    sigma_y = sigma_y+sigma_u;
143end
144
145fact = norminv((1-options_.forecasts.conf_sig)/2,0,1);
146if nargout==3
147    var_yf_ME=var_yf;
148    var_yf_ME(:,options_.varobs_id)=var_yf(:,options_.varobs_id)+repmat(diag(M_.H)',horizon,1);
149    int_width_ME = zeros(horizon,M_.endo_nbr);
150end
151
152int_width = zeros(iter,nvar);
153for i=1:nvar
154    int_width(:,i) = fact*sqrt(var_yf(:,i));
155    if nargout==3
156        int_width_ME(:,i) = -fact*sqrt(var_yf_ME(:,i));
157    end
158end
159