stepan
/
dynare
1
0
Fork 0
Code Issues Pull Requests Packages Projects Releases Wiki Activity

prior_posterior_statistics_core.m: consolidate both forecast functions into one inline function

kalman-mex
Johannes Pfeifer 2023-09-16 13:44:26 +02:00 committed by Sébastien Villemot
parent dc9560e610
commit bf7ac27fd7
No known key found for this signature in database
GPG Key ID: 2CECE9350ECEBE4A
3 changed files with 55 additions and 119 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

67
matlab/forcst2.m
View File

@ -1,67 +0,0 @@
function yf=forcst2(y0,horizon,dr,n)
% function yf=forcst2(y0,horizon,dr,n)
%
% computes forecasts based on first order model solution, given shocks
% drawn from the shock distribution, but not including measurement error
% Inputs:
% - y0 [endo_nbr by maximum_endo_lag] matrix of starting values
% - dr [structure] structure with Dynare decision rules
% - e [horizon by exo_nbr] matrix with shocks
% - n [scalar] number of repetitions
%
% Outputs:
% - yf [horizon+ykmin_ by endo_nbr by n] array of forecasts
%
% Copyright © 2008-2017 Dynare Team
%
% 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
% along with Dynare. If not, see <https://www.gnu.org/licenses/>.
global M_ options_
Sigma_e_ = M_.Sigma_e;
endo_nbr = M_.endo_nbr;
exo_nbr = M_.exo_nbr;
ykmin_ = M_.maximum_endo_lag;
k1 = [ykmin_:-1:1];
k2 = dr.kstate(find(dr.kstate(:,2) <= ykmin_+1),[1 2]);
k2 = k2(:,1)+(ykmin_+1-k2(:,2))*endo_nbr;
% eliminate shocks with 0 variance
i_exo_var = setdiff([1:exo_nbr],find(diag(Sigma_e_) == 0));
nxs = length(i_exo_var);
chol_S = chol(Sigma_e_(i_exo_var,i_exo_var));
if ~isempty(Sigma_e_)
e = randn(nxs,n,horizon);
end
B1 = dr.ghu(:,i_exo_var)*chol_S';
yf = zeros(endo_nbr,horizon+ykmin_,n);
yf(:,1:ykmin_,:,:) = repmat(y0,[1,1,n]);
j = ykmin_*endo_nbr;
for i=ykmin_+(1:horizon)
tempx1 = reshape(yf(:,k1,:),[j,n]);
tempx = tempx1(k2,:);
yf(:,i,:) = dr.ghx*tempx+B1*squeeze(e(:,:,i-ykmin_));
k1 = k1+1;
end
yf(dr.order_var,:,:) = yf;
yf=permute(yf,[2 1 3]);

50
matlab/forcst2a.m
View File

@ -1,50 +0,0 @@
function yf=forcst2a(y0,dr,e)
% function yf=forcst2a(y0,dr,e)
% computes forecasts based on first order model solution, assuming the absence of shocks
% Inputs:
% - y0 [endo_nbr by maximum_endo_lag] matrix of starting values
% - dr [structure] structure with Dynare decision rules
% - e [horizon by exo_nbr] matrix with shocks
%
% Outputs:
% - yf [horizon+maximum_endo_lag,endo_nbr] matrix of forecasts
%
% Copyright © 2008-2017 Dynare Team
%
% 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
% along with Dynare. If not, see <https://www.gnu.org/licenses/>.
global M_ options_
endo_nbr = M_.endo_nbr;
ykmin_ = M_.maximum_endo_lag;
horizon = size(e,1);
k1 = [ykmin_:-1:1];
k2 = dr.kstate(find(dr.kstate(:,2) <= ykmin_+1),[1 2]);
k2 = k2(:,1)+(ykmin_+1-k2(:,2))*endo_nbr;
yf = zeros(horizon+ykmin_,endo_nbr);
yf(1:ykmin_,:) = y0';
j = ykmin_*endo_nbr;
for i=ykmin_+(1:horizon)
tempx = yf(k1,:)';
yf(i,:) = tempx(k2)'*dr.ghx';
k1 = k1+1;
end
yf(:,dr.order_var) = yf;

57
matlab/prior_posterior_statistics_core.m
View File

@ -314,7 +314,7 @@ for b=fpar:B
end
if horizon
yyyy = alphahat(iendo,i_last_obs);
yf = forcst2a(yyyy,dr,zeros(horizon,exo_nbr));
yf = simulate_posterior_forecasts(yyyy,dr,horizon,false,M_.Sigma_e,1);
if options_.prefilter
% add mean
yf(:,IdObs) = yf(:,IdObs)+repmat(mean_varobs, ...
@ -331,7 +331,7 @@ for b=fpar:B
else
yf = yf+repmat(SteadyState',horizon+maxlag,1);
end
yf1 = forcst2(yyyy,horizon,dr,1);
yf1 = simulate_posterior_forecasts(yyyy,dr,horizon,false,M_.Sigma_e,1);
if options_.prefilter == 1
% add mean
yf1(:,IdObs,:) = yf1(:,IdObs,:)+ ...
@ -543,3 +543,56 @@ if RemoteFlag==1
end
dyn_waitbar_close(h);
function yf=simulate_posterior_forecasts(y0,dr,horizon,stochastic_indicator,Sigma_e,n)
% function yf=forcst2(y0,horizon,dr,n)
%
% computes forecasts based on first order model solution, given shocks
% drawn from the shock distribution, but not including measurement error
% Inputs:
% - y0 [endo_nbr by maximum_endo_lag] matrix of starting values
% - dr [structure] structure with Dynare decision rules
% - horizon [scalar] number of forecast periods
% - stochastic_indicator [boolean] indicator whether to consider stochastic shocks
% - Sigma_e [integer] covariance matrix of shocks
% - n [scalar] number of repetitions
%
% Outputs:
% - yf [horizon+ykmin_ by endo_nbr by n] array of forecasts
if nargin< 4
stochastic_indicator=false;
n=1;
end
%select states
k2 = dr.inv_order_var(dr.state_var);
if stochastic_indicator
% eliminate shocks with 0 variance
i_exo_var = setdiff(1:length(Sigma_e),find(diag(Sigma_e) == 0));
nxs = length(i_exo_var);
chol_S = chol(Sigma_e(i_exo_var,i_exo_var));
if ~isempty(Sigma_e)
e = randn(nxs,n,horizon);
end
B1 = dr.ghu(:,i_exo_var)*chol_S';
end
endo_nbr=length(y0);
yf = zeros(endo_nbr,1+horizon,n);
yf(:,1,:,:) = repmat(y0,[1,1,n]);
for iter=1:horizon
if stochastic_indicator
yf(:,iter+1,:) = dr.ghx*squeeze(yf(k2,iter,:))+B1*squeeze(e(:,:,iter));
else
yf(:,iter+1,:) = dr.ghx*squeeze(yf(k2,iter,:));
end
end
yf(dr.order_var,:,:) = yf;
yf=permute(yf,[2 1 3]);