2016-06-19 12:41:51 +02:00
|
|
|
function [y_,int_width,int_width_ME]=simultxdet(y0,ex,ex_det, iorder,var_list,M_,oo_,options_)
|
2009-10-29 18:16:10 +01:00
|
|
|
%function [y_,int_width]=simultxdet(y0,ex,ex_det, iorder,var_list,M_,oo_,options_)
|
|
|
|
%
|
|
|
|
% Simulates a stochastic model in the presence of deterministic exogenous shocks
|
|
|
|
%
|
|
|
|
% INPUTS:
|
|
|
|
% y0: initial values, of length M_.maximum_lag
|
|
|
|
% ex: matrix of stochastic exogenous shocks, starting at period 1
|
|
|
|
% ex_det: matrix of deterministic exogenous shocks, starting at period 1-M_.maximum_lag
|
|
|
|
% iorder: order of approximation
|
|
|
|
% var_list: list of endogenous variables to simulate
|
2016-06-19 12:41:51 +02:00
|
|
|
% int_width_ME:distance between upper bound and
|
|
|
|
% mean forecast when considering measurement error
|
|
|
|
% OUTPUTS:
|
|
|
|
% yf: mean forecast
|
|
|
|
% int_width: distance between upper bound and
|
|
|
|
% mean forecast
|
|
|
|
% int_width_ME:distance between upper bound and
|
|
|
|
% mean forecast when considering measurement error
|
2009-10-29 18:16:10 +01:00
|
|
|
%
|
|
|
|
% The forecast horizon is equal to size(ex, 1).
|
|
|
|
% The condition size(ex,1)+M_.maximum_lag=size(ex_det,1) must be verified
|
|
|
|
% for consistency.
|
|
|
|
|
2017-10-10 10:05:59 +02:00
|
|
|
% Copyright (C) 2008-2018 Dynare Team
|
2009-10-29 18:16:10 +01:00
|
|
|
%
|
|
|
|
% This file is part of Dynare.
|
|
|
|
%
|
|
|
|
% Dynare is free software: you can redistribute it and/or modify
|
|
|
|
% it under the terms of the GNU General Public License as published by
|
|
|
|
% the Free Software Foundation, either version 3 of the License, or
|
|
|
|
% (at your option) any later version.
|
|
|
|
%
|
|
|
|
% Dynare is distributed in the hope that it will be useful,
|
|
|
|
% but WITHOUT ANY WARRANTY; without even the implied warranty of
|
|
|
|
% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
|
|
|
|
% GNU General Public License for more details.
|
|
|
|
%
|
|
|
|
% You should have received a copy of the GNU General Public License
|
2021-06-09 17:33:48 +02:00
|
|
|
% along with Dynare. If not, see <https://www.gnu.org/licenses/>.
|
2009-10-29 18:16:10 +01:00
|
|
|
|
2009-12-16 18:17:34 +01:00
|
|
|
dr = oo_.dr;
|
|
|
|
ykmin = M_.maximum_lag;
|
|
|
|
endo_nbr = M_.endo_nbr;
|
|
|
|
exo_det_steady_state = oo_.exo_det_steady_state;
|
2012-11-16 20:05:13 +01:00
|
|
|
nstatic = M_.nstatic;
|
|
|
|
nspred = M_.nspred;
|
2009-12-16 18:17:34 +01:00
|
|
|
nc = size(dr.ghx,2);
|
|
|
|
iter = size(ex,1);
|
|
|
|
if size(ex_det, 1) ~= iter+ykmin
|
|
|
|
error('Size mismatch: number of forecasting periods for stochastic exogenous and deterministic exogenous don''t match')
|
|
|
|
end
|
|
|
|
nx = size(dr.ghu,2);
|
|
|
|
y_ = zeros(size(y0,1),iter+ykmin);
|
|
|
|
y_(:,1:ykmin) = y0;
|
|
|
|
k1 = [ykmin:-1:1];
|
|
|
|
k2 = dr.kstate(find(dr.kstate(:,2) <= ykmin+1),[1 2]);
|
|
|
|
k2 = k2(:,1)+(ykmin+1-k2(:,2))*endo_nbr;
|
|
|
|
k3 = M_.lead_lag_incidence(1:ykmin,:)';
|
|
|
|
k3 = find(k3(:));
|
|
|
|
k4 = dr.kstate(find(dr.kstate(:,2) < ykmin+1),[1 2]);
|
|
|
|
k4 = k4(:,1)+(ykmin+1-k4(:,2))*endo_nbr;
|
|
|
|
|
2017-10-10 10:05:59 +02:00
|
|
|
nvar = length(var_list);
|
2009-12-16 18:17:34 +01:00
|
|
|
if nvar == 0
|
2009-10-29 18:16:10 +01:00
|
|
|
nvar = endo_nbr;
|
|
|
|
ivar = [1:nvar];
|
2009-12-16 18:17:34 +01:00
|
|
|
else
|
2009-10-29 18:16:10 +01:00
|
|
|
ivar=zeros(nvar,1);
|
|
|
|
for i=1:nvar
|
2017-10-10 10:05:59 +02:00
|
|
|
i_tmp = strmatch(var_list{i}, M_.endo_names, 'exact');
|
2009-12-16 18:17:34 +01:00
|
|
|
if isempty(i_tmp)
|
2017-10-10 10:05:59 +02:00
|
|
|
disp(var_list{i})
|
2009-12-16 18:17:34 +01:00
|
|
|
error (['One of the variable specified does not exist']) ;
|
|
|
|
else
|
|
|
|
ivar(i) = i_tmp;
|
|
|
|
end
|
2009-10-29 18:16:10 +01:00
|
|
|
end
|
2009-12-16 18:17:34 +01:00
|
|
|
end
|
2009-10-29 18:16:10 +01:00
|
|
|
|
2009-12-16 18:17:34 +01:00
|
|
|
if iorder == 1
|
2009-10-29 18:16:10 +01:00
|
|
|
for i = ykmin+1: iter+ykmin
|
2009-12-16 18:17:34 +01:00
|
|
|
tempx1 = y_(dr.order_var,k1);
|
|
|
|
tempx2 = tempx1-repmat(dr.ys(dr.order_var),1,ykmin);
|
|
|
|
tempx = tempx2(k2);
|
2017-10-10 10:05:59 +02:00
|
|
|
y_(dr.order_var,i) = dr.ys(dr.order_var)+dr.ghx*tempx+dr.ghu*ex(i-ykmin,:)';
|
2010-01-05 11:46:10 +01:00
|
|
|
for j=1:min(ykmin+M_.exo_det_length+1-i,M_.exo_det_length)
|
2009-10-29 18:16:10 +01:00
|
|
|
y_(dr.order_var,i) = y_(dr.order_var,i) + dr.ghud{j}*(ex_det(i+j-1,:)'-exo_det_steady_state);
|
2010-01-05 11:46:10 +01:00
|
|
|
end
|
2017-05-16 15:10:20 +02:00
|
|
|
|
2010-01-05 11:46:10 +01:00
|
|
|
k1 = k1+1;
|
2009-10-29 18:16:10 +01:00
|
|
|
end
|
2009-12-16 18:17:34 +01:00
|
|
|
elseif iorder == 2
|
2009-10-29 18:16:10 +01:00
|
|
|
for i = ykmin+1: iter+ykmin
|
2009-12-16 18:17:34 +01:00
|
|
|
tempx1 = y_(dr.order_var,k1);
|
|
|
|
tempx2 = tempx1-repmat(dr.ys(dr.order_var),1,ykmin);
|
|
|
|
tempx = tempx2(k2);
|
|
|
|
tempu = ex(i-ykmin,:)';
|
|
|
|
tempuu = kron(tempu,tempu);
|
2010-01-05 11:46:10 +01:00
|
|
|
tempxx = kron(tempx,tempx);
|
|
|
|
tempxu = kron(tempx,tempu);
|
|
|
|
y_(dr.order_var,i) = dr.ys(dr.order_var)+dr.ghs2/2+dr.ghx*tempx+ ...
|
|
|
|
dr.ghu*tempu+0.5*(dr.ghxx*tempxx+dr.ghuu*tempuu)+dr.ghxu* ...
|
|
|
|
tempxu;
|
|
|
|
for j=1:min(ykmin+M_.exo_det_length+1-i,M_.exo_det_length)
|
2009-12-16 18:17:34 +01:00
|
|
|
tempud = ex_det(i+j-1,:)'-exo_det_steady_state;
|
|
|
|
tempudud = kron(tempud,tempud);
|
|
|
|
tempxud = kron(tempx,tempud);
|
|
|
|
tempuud = kron(tempu,tempud);
|
|
|
|
y_(dr.order_var,i) = y_(dr.order_var,i) + dr.ghud{j}*tempud + ...
|
|
|
|
dr.ghxud{j}*tempxud + dr.ghuud{j}*tempuud + ...
|
|
|
|
0.5*dr.ghudud{j,j}*tempudud;
|
|
|
|
for k=1:j-1
|
|
|
|
tempudk = ex_det(i+k-1,:)'-exo_det_steady_state;
|
|
|
|
tempududk = kron(tempudk,tempud);
|
|
|
|
y_(dr.order_var,i) = y_(dr.order_var,i) + ...
|
|
|
|
dr.ghudud{k,j}*tempududk;
|
|
|
|
end
|
2010-01-05 11:46:10 +01:00
|
|
|
end
|
|
|
|
k1 = k1+1;
|
2009-10-29 18:16:10 +01:00
|
|
|
end
|
2009-12-16 18:17:34 +01:00
|
|
|
end
|
|
|
|
|
2012-11-16 20:05:13 +01:00
|
|
|
[A,B] = kalman_transition_matrix(dr,nstatic+(1:nspred),1:nc,M_.exo_nbr);
|
2009-10-29 18:16:10 +01:00
|
|
|
|
2009-12-16 18:17:34 +01:00
|
|
|
inv_order_var = dr.inv_order_var;
|
|
|
|
ghx1 = dr.ghx(inv_order_var(ivar),:);
|
|
|
|
ghu1 = dr.ghu(inv_order_var(ivar),:);
|
2009-10-29 18:16:10 +01:00
|
|
|
|
2009-12-16 18:17:34 +01:00
|
|
|
sigma_u = B*M_.Sigma_e*B';
|
|
|
|
sigma_u1 = ghu1*M_.Sigma_e*ghu1';
|
|
|
|
sigma_y = 0;
|
2009-10-29 18:16:10 +01:00
|
|
|
|
2016-06-19 14:31:19 +02:00
|
|
|
var_yf=NaN(iter,nvar); %initialize
|
2009-12-16 18:17:34 +01:00
|
|
|
for i=1:iter
|
|
|
|
sigma_y1 = ghx1*sigma_y*ghx1'+sigma_u1;
|
|
|
|
var_yf(i,:) = diag(sigma_y1)';
|
|
|
|
if i == iter
|
|
|
|
break
|
|
|
|
end
|
|
|
|
sigma_u = A*sigma_u*A';
|
|
|
|
sigma_y = sigma_y+sigma_u;
|
|
|
|
end
|
2009-10-29 18:16:10 +01:00
|
|
|
|
2016-07-21 12:17:30 +02:00
|
|
|
fact = norminv((1-options_.forecasts.conf_sig)/2,0,1);
|
2016-06-19 12:41:51 +02:00
|
|
|
if nargout==3
|
|
|
|
var_yf_ME=var_yf;
|
|
|
|
var_yf_ME(:,options_.varobs_id)=var_yf(:,options_.varobs_id)+repmat(diag(M_.H)',horizon,1);
|
|
|
|
int_width_ME = zeros(horizon,M_.endo_nbr);
|
|
|
|
end
|
2009-10-29 18:16:10 +01:00
|
|
|
|
2016-06-19 12:41:51 +02:00
|
|
|
int_width = zeros(iter,nvar);
|
2009-12-16 18:17:34 +01:00
|
|
|
for i=1:nvar
|
2009-10-29 18:16:10 +01:00
|
|
|
int_width(:,i) = fact*sqrt(var_yf(:,i));
|
2017-05-16 15:10:20 +02:00
|
|
|
if nargout==3
|
2016-06-19 12:41:51 +02:00
|
|
|
int_width_ME(:,i) = -fact*sqrt(var_yf_ME(:,i));
|
2017-05-16 15:10:20 +02:00
|
|
|
end
|
2009-12-16 18:17:34 +01:00
|
|
|
end
|