v4: added forecast out of posterior mode (or maximum likelihood, non-tested)
git-svn-id: https://www.dynare.org/svn/dynare/dynare_v4@1710 ac1d8469-bf42-47a9-8791-bf33cf982152time-shift
michel
parent
5572943251
commit
bcf9da214e
|
|
@ -1389,7 +1389,7 @@ if (~((any(bayestopt_.pshape > 0) & options_.mh_replic) | (any(bayestopt_.pshape
|
|||
end
|
||||
|
||||
if options_.forecast > 0 & options_.mh_replic == 0 & ~options_.load_mh_file
|
||||
forecast(var_list);
|
||||
forecast(var_list_,'smoother');
|
||||
end
|
||||
|
||||
pindx = estim_params_.param_vals(:,1);
|
||||
|
|
|
|||
125
matlab/forcst.m
125
matlab/forcst.m
|
|
@ -1,54 +1,77 @@
|
|||
function [yf,int_width]=forcst(dr,y0,k,var_list)
|
||||
global M_ oo_ options_
|
||||
|
||||
make_ex_;
|
||||
yf = simult_(y0,dr,zeros(k,M_.exo_nbr),1);
|
||||
nstatic = dr.nstatic;
|
||||
npred = dr.npred;
|
||||
nc = size(dr.ghx,2);
|
||||
endo_nbr = M_.endo_nbr;
|
||||
inv_order_var = dr.inv_order_var;
|
||||
[A,B] = kalman_transition_matrix(dr,nstatic+(1:npred),1:nc,dr.transition_auxiliary_variables);
|
||||
|
||||
nvar = size(var_list,1);
|
||||
if nvar == 0
|
||||
nvar = M_.endo_nbr;
|
||||
ivar = [1:nvar];
|
||||
else
|
||||
ivar=zeros(nvar,1);
|
||||
function [yf,int_width]=forcst(dr,y0,horizon,var_list)
|
||||
|
||||
% function [yf,int_width]=forecst(dr,y0,horizon,var_list)
|
||||
% computes mean forecast for a given value of the parameters
|
||||
% computes also confidence band for the forecast
|
||||
%
|
||||
% INPUTS:
|
||||
% dr: structure containing decision rules
|
||||
% y0: initial values
|
||||
% horizon: nbr of periods to forecast
|
||||
% var_list: list of variables (character matrix)
|
||||
%
|
||||
% OUTPUTS:
|
||||
% yf: mean forecast
|
||||
% int_width: distance between upper bound and
|
||||
% mean forecast
|
||||
%
|
||||
% SPECIAL REQUIREMENTS
|
||||
% none
|
||||
%
|
||||
% part of DYNARE, copyright Dynare Team (2003-2008)
|
||||
% Gnu Public License.
|
||||
|
||||
global M_ oo_ options_
|
||||
|
||||
make_ex_;
|
||||
yf = simult_(y0,dr,zeros(horizon,M_.exo_nbr),1);
|
||||
nstatic = dr.nstatic;
|
||||
npred = dr.npred;
|
||||
nc = size(dr.ghx,2);
|
||||
endo_nbr = M_.endo_nbr;
|
||||
inv_order_var = dr.inv_order_var;
|
||||
[A,B] = kalman_transition_matrix(dr,nstatic+(1:npred),1:nc,dr.transition_auxiliary_variables);
|
||||
|
||||
nvar = size(var_list,1);
|
||||
if nvar == 0
|
||||
nvar = M_.endo_nbr;
|
||||
ivar = [1:nvar];
|
||||
else
|
||||
ivar=zeros(nvar,1);
|
||||
for i=1:nvar
|
||||
i_tmp = strmatch(var_list(i,:),M_.endo_names,'exact');
|
||||
if isempty(i_tmp)
|
||||
disp(var_list(i,:));
|
||||
error (['One of the variable specified does not exist']) ;
|
||||
else
|
||||
ivar(i) = i_tmp;
|
||||
end
|
||||
end
|
||||
end
|
||||
|
||||
ghx1 = dr.ghx(inv_order_var(ivar),:);
|
||||
ghu1 = dr.ghu(inv_order_var(ivar),:);
|
||||
|
||||
sigma_u = B*M_.Sigma_e*B';
|
||||
sigma_u1 = ghu1*M_.Sigma_e*ghu1';
|
||||
sigma_y = 0;
|
||||
|
||||
for i=1:horizon
|
||||
sigma_y1 = ghx1*sigma_y*ghx1'+sigma_u1;
|
||||
var_yf(i,:) = diag(sigma_y1)';
|
||||
if i == horizon
|
||||
break
|
||||
end
|
||||
sigma_u = A*sigma_u*A';
|
||||
sigma_y = sigma_y+sigma_u;
|
||||
end
|
||||
|
||||
fact = qnorm((1-options_.conf_sig)/2,0,1);
|
||||
|
||||
int_width = zeros(horizon,M_.endo_nbr);
|
||||
for i=1:nvar
|
||||
i_tmp = strmatch(var_list(i,:),M_.endo_names,'exact');
|
||||
if isempty(i_tmp)
|
||||
disp(var_list(i,:));
|
||||
error (['One of the variable specified does not exist']) ;
|
||||
else
|
||||
ivar(i) = i_tmp;
|
||||
end
|
||||
int_width(:,i) = fact*sqrt(var_yf(:,i));
|
||||
end
|
||||
end
|
||||
|
||||
ghx1 = dr.ghx(inv_order_var(ivar),:);
|
||||
ghu1 = dr.ghu(inv_order_var(ivar),:);
|
||||
|
||||
sigma_u = B*M_.Sigma_e*B';
|
||||
sigma_u1 = ghu1*M_.Sigma_e*ghu1';
|
||||
sigma_y = 0;
|
||||
|
||||
for i=1:k
|
||||
sigma_y1 = ghx1*sigma_y*ghx1'+sigma_u1;
|
||||
var_yf(i,:) = diag(sigma_y1)';
|
||||
if i == k
|
||||
break
|
||||
end
|
||||
sigma_u = A*sigma_u*A';
|
||||
sigma_y = sigma_y+sigma_u;
|
||||
end
|
||||
|
||||
fact = qnorm((1-options_.conf_sig)/2,0,1);
|
||||
|
||||
int_width = zeros(k,M_.endo_nbr);
|
||||
for i=1:nvar
|
||||
int_width(:,i) = fact*sqrt(var_yf(:,i));
|
||||
end
|
||||
|
||||
yf = yf(ivar,:);
|
||||
yf = yf(ivar,:);
|
||||
|
||||
|
|
@ -1,49 +1,125 @@
|
|||
% Copyright (C) 2005 Michel Juillard
|
||||
%
|
||||
function forecast(var_list)
|
||||
global options_ dr_ oo_ M_
|
||||
|
||||
old_options = options_;
|
||||
options_ = set_default_option(options_,'periods',40);
|
||||
if options_.periods == 0
|
||||
options_.periods = 40;
|
||||
end
|
||||
options_ = set_default_option(options_,'conf_sig',0.9);
|
||||
|
||||
if size(oo_.endo_simul,2) < M_.maximum_lag
|
||||
y0 = repmat(oo_.steady_state,M_.maximum_lag);
|
||||
else
|
||||
y0 = oo_.endo_simul(:,1:M_.maximum_lag);
|
||||
end
|
||||
|
||||
if M_.exo_det_nbr == 0
|
||||
[yf,int_width] = forcst(oo_.dr,y0,options_.periods,var_list);
|
||||
else
|
||||
exo_det_length = size(oo_.exo_det_simul,1);
|
||||
if options_.periods > exo_det_length
|
||||
ex = zeros(options_.periods,M_.exo_nbr);
|
||||
oo_.exo_det_simul = [ oo_.exo_det_simul;...
|
||||
repmat(oo_.exo_det_steady_state',...
|
||||
options_.periods- ...
|
||||
exo_det_length,1)];
|
||||
%ex_det_length,1),1)];
|
||||
elseif options_.periods < exo_det_length
|
||||
ex = zeros(exo_det_length,M_.exo_nbr);
|
||||
end
|
||||
[yf,int_width] = simultxdet(y0,dr_,ex,oo_.exo_det_simul,...
|
||||
options_.order,var_list);
|
||||
end
|
||||
|
||||
for i=1:M_.endo_nbr
|
||||
eval(['oo_.forecast.Mean.' M_.endo_names(i,:) '= yf(' int2str(i) ',M_.maximum_lag+1:end)'';']);
|
||||
eval(['oo_.forecast.HPDinf.' M_.endo_names(i,:) '= yf(' int2str(i) ',M_.maximum_lag+1:end)''-' ...
|
||||
' int_width(:,' int2str(i) ');']);
|
||||
eval(['oo_.forecast.HPDsup.' M_.endo_names(i,:) '= yf(' int2str(i) ',M_.maximum_lag+1:end)''+' ...
|
||||
' int_width(:,' int2str(i) ');']);
|
||||
end
|
||||
function forecast(var_list,task)
|
||||
|
||||
for i=1:M_.exo_det_nbr
|
||||
eval(['oo_.forecast.Exogenous.' lgx_det_(i,:) '= M_.exo_det_simul(:,' int2str(i) ');']);
|
||||
end
|
||||
|
||||
options_ = old_options;
|
||||
% function forecast(var_list,task)
|
||||
% computes mean forecast for a given value of the parameters
|
||||
% computes also confidence band for the forecast
|
||||
%
|
||||
% INPUTS
|
||||
% var_list: list of variables (character matrix)
|
||||
% task: indicates how to initialize the forecast
|
||||
% either 'simul' or 'smoother'
|
||||
% OUTPUTS
|
||||
% nothing is returned but the procedure saves output
|
||||
% in oo_.forecast.Mean
|
||||
% oo_.forecast.HPDinf
|
||||
% oo_.forecast.HPDsup
|
||||
%
|
||||
% SPECIAL REQUIREMENTS
|
||||
% none
|
||||
%
|
||||
% part of DYNARE, copyright Dynare Team (2003-2008)
|
||||
% Gnu Public License.
|
||||
|
||||
global options_ dr_ oo_ M_
|
||||
|
||||
old_options = options_;
|
||||
|
||||
maximum_lag = M_.maximum_lag;
|
||||
y_smoothed = oo_.SmoothedVariables;
|
||||
options_ = set_default_option(options_,'periods',40);
|
||||
if options_.periods == 0
|
||||
options_.periods = 40;
|
||||
end
|
||||
horizon = options_.periods;
|
||||
options_ = set_default_option(options_,'conf_sig',0.9);
|
||||
|
||||
endo_names = M_.endo_names;
|
||||
if isempty(var_list)
|
||||
varlist = endo_names;
|
||||
i_var = 1:M_.endo_nbr;
|
||||
else
|
||||
i_var = [];
|
||||
for i = 1:size(var_list)
|
||||
tmp = strmatch(var_list(i,:),endo_names,'exact');
|
||||
if isempty(tmp)
|
||||
error([var_list(i,:) ' isn''t and endogenous variable'])
|
||||
end
|
||||
i_var = [i_var; tmp];
|
||||
end
|
||||
end
|
||||
n_var = length(i_var);
|
||||
|
||||
trend = 0;
|
||||
switch task
|
||||
case 'simul'
|
||||
if size(oo_.endo_simul,2) < maximum_lag
|
||||
y0 = repmat(oo_.steady_state,1,maximum_lag);
|
||||
else
|
||||
y0 = oo_.endo_simul(:,1:maximum_lag);
|
||||
end
|
||||
case 'smoother'
|
||||
y0 = zeros(M_.endo_nbr,maximum_lag);
|
||||
for i = 1:M_.endo_nbr
|
||||
v_name = deblank(M_.endo_names(i,:));
|
||||
y0(i,:) = y_smoothed.(v_name)(end-maximum_lag+1:end)+oo_.dr.ys(i);
|
||||
end
|
||||
gend = size(y0,2);
|
||||
if isfield(oo_.Smoother,'TrendCoeffs')
|
||||
var_obs = options_.varobs;
|
||||
endo_names = M_.endo_names;
|
||||
order_var = oo_.dr.order_var;
|
||||
i_var_obs = [];
|
||||
trend_coeffs = [];
|
||||
for i=1:size(var_obs,1)
|
||||
tmp = strmatch(var_obs(i,:),endo_names(i_var,:),'exact');
|
||||
if ~isempty(tmp)
|
||||
i_var_obs = [ i_var_obs; tmp];
|
||||
trend_coeffs = [trend_coeffs; oo_.Smoother.TrendCoeffs(i)];
|
||||
end
|
||||
end
|
||||
trend = trend_coeffs*(gend+(1-M_.maximum_lag:horizon));
|
||||
end
|
||||
otherwise
|
||||
error('Wrong flag value')
|
||||
end
|
||||
|
||||
if M_.exo_det_nbr == 0
|
||||
[yf,int_width] = forcst(oo_.dr,y0,options_.periods,var_list);
|
||||
else
|
||||
exo_det_length = size(oo_.exo_det_simul,1);
|
||||
if options_.periods > exo_det_length
|
||||
ex = zeros(options_.periods,M_.exo_nbr);
|
||||
oo_.exo_det_simul = [ oo_.exo_det_simul;...
|
||||
repmat(oo_.exo_det_steady_state',...
|
||||
options_.periods- ...
|
||||
exo_det_length,1)];
|
||||
%ex_det_length,1),1)];
|
||||
elseif options_.periods < exo_det_length
|
||||
ex = zeros(exo_det_length,M_.exo_nbr);
|
||||
end
|
||||
[yf,int_width] = simultxdet(y0,dr_,ex,oo_.exo_det_simul,...
|
||||
options_.order,var_list);
|
||||
end
|
||||
|
||||
if ~isscalar(trend)
|
||||
yf(i_var_obs,:) = yf(i_var_obs,:) + trend;
|
||||
end
|
||||
|
||||
for i=1:n_var
|
||||
eval(['oo_.forecast.Mean.' var_list(i,:) '= yf(' int2str(i) ',maximum_lag+1:end)'';']);
|
||||
eval(['oo_.forecast.HPDinf.' var_list(i,:) '= yf(' int2str(i) ',maximum_lag+1:end)''-' ...
|
||||
' int_width(:,' int2str(i) ');']);
|
||||
eval(['oo_.forecast.HPDsup.' var_list(i,:) '= yf(' int2str(i) ',maximum_lag+1:end)''+' ...
|
||||
' int_width(:,' int2str(i) ');']);
|
||||
end
|
||||
|
||||
for i=1:M_.exo_det_nbr
|
||||
eval(['oo_.forecast.Exogenous.' lgx_det_(i,:) '= M_.exo_det_simul(:,' int2str(i) ');']);
|
||||
end
|
||||
|
||||
if options_.nograph == 0
|
||||
forecast_graphs(var_list);
|
||||
end
|
||||
|
||||
options_ = old_options;
|
||||
|
||||
|
|
|
|||
|
|
@ -0,0 +1,77 @@
|
|||
function forecast_graphs(var_list)
|
||||
|
||||
global M_ oo_ options_
|
||||
|
||||
nc = 4;
|
||||
nr = 3;
|
||||
exo_nbr = M_.exo_nbr;
|
||||
endo_names = M_.endo_names;
|
||||
fname = M_.fname;
|
||||
% $$$ varobs = options_.varobs;
|
||||
% $$$ y = oo_.SmoothedVariables;
|
||||
% $$$ ys = oo_.dr.ys;
|
||||
% $$$ gend = size(y,2);
|
||||
yf = oo_.forecast.Mean;
|
||||
hpdinf = oo_.forecast.HPDinf;
|
||||
hpdsup = oo_.forecast.HPDsup;
|
||||
if isempty(var_list)
|
||||
varlist = endo_names;
|
||||
i_var = 1:M_.endo_nbr;
|
||||
else
|
||||
i_var = [];
|
||||
for i = 1:size(var_list)
|
||||
tmp = strmatch(var_list(i,:),endo_names,'exact');
|
||||
if isempty(tmp)
|
||||
error([var_list(i,:) ' isn''t and endogenous variable'])
|
||||
end
|
||||
i_var = [i_var; tmp];
|
||||
end
|
||||
end
|
||||
nvar = length(i_var);
|
||||
|
||||
% $$$ % build trend for smoothed variables if necessary
|
||||
% $$$ trend = zeros(size(varobs,1),10);
|
||||
% $$$ if isfield(oo_.Smoother,'TrendCoeffs')
|
||||
% $$$ trend_coeffs = oo_.Smoother.TrendCoeffs;
|
||||
% $$$ trend = trend_coeffs*(gend-9:gend);
|
||||
% $$$ end
|
||||
|
||||
% create subdirectory <fname>/graphs if id doesn't exist
|
||||
if ~exist([fname '/graphs'])
|
||||
mkdir(fname,'graphs');
|
||||
end
|
||||
|
||||
m = 1;
|
||||
n_fig = 1;
|
||||
figure('Name','Forecasts (I)')
|
||||
for j= 1:nvar
|
||||
if m > nc*nr;
|
||||
eval(['print -depsc ' fname '/graphs/forcst' int2str(n_fig) ...
|
||||
'.eps;'])
|
||||
n_fig =n_fig+1;
|
||||
eval(['figure(''Name'',''Forecast (' int2str(n_fig) ')'');']);
|
||||
m = 1;
|
||||
end
|
||||
subplot(nr,nc,m);
|
||||
vn = deblank(endo_names(i_var(j),:));
|
||||
obs = 0;
|
||||
% $$$ k = strmatch(vn,varobs,'exact');
|
||||
% $$$ if ~isempty(k)
|
||||
% $$$ yy = y.(vn)(end-9:end) + repmat(ys(i_var(j)),10,1)+trend(k,:)';
|
||||
% $$$ plot(yy);
|
||||
% $$$ hold on
|
||||
% $$$ obs = 10;
|
||||
% $$$ end
|
||||
plot([NaN(obs,1); yf.(vn)]);
|
||||
hold on
|
||||
plot([NaN(obs,1); hpdinf.(vn)]);
|
||||
hold on
|
||||
plot([NaN(obs,1); hpdsup.(vn)]);
|
||||
title(vn,'Interpreter','none');
|
||||
hold off
|
||||
m = m + 1;
|
||||
end
|
||||
|
||||
if m > 1
|
||||
eval(['print -deps ' fname '/graphs/forcst' int2str(n_fig) '.eps;'])
|
||||
end
|
||||
Loading…
Reference in New Issue