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