diff --git a/matlab/imcforecast.m b/matlab/imcforecast.m index b3355d7ec..fba682a51 100644 --- a/matlab/imcforecast.m +++ b/matlab/imcforecast.m @@ -1,6 +1,32 @@ -function imcforecast(ptype,cV,cS,cL,H,mcValue,B,ci) +function imcforecast(constrained_paths, constrained_vars, options_cond_fcst) +% Computes conditional forecasts. +% +% INPUTS +% o constrained_paths [double] m*p array, where m is the number of constrained endogenous variables and p is the number of constrained periods. +% o constrained_vars [char] m*x array holding the names of the controlled endogenous variables. +% o options_cond_fcst [structure] containing the options. The fields are: +% + replic [integer] scalar, number of monte carlo simulations. +% + parameter_set [char] values of the estimated parameters: +% "posterior_mode", +% "posterior_mean", +% "posterior_median", +% "prior_mode" or +% "prior mean". +% [double] np*1 array, values of the estimated parameters. +% + controlled_varexo [char] m*x array, list of controlled exogenous variables. +% + conf_sig [double] scalar in [0,1], probability mass covered by the confidence bands. +% +% OUTPUTS +% None. +% +% SPECIAL REQUIREMENTS +% This routine has to be called after an estimation statement or an estimated_params block. +% +% REMARKS +% [1] Results are stored in a structure which is saved in a mat file called conditional_forecasts.mat. +% [2] Use the function plot_icforecast to plot the results. -% Copyright (C) 2006 Dynare Team +% Copyright (C) 2006-2009 Dynare Team % % This file is part of Dynare. % @@ -17,33 +43,97 @@ function imcforecast(ptype,cV,cS,cL,H,mcValue,B,ci) % You should have received a copy of the GNU General Public License % along with Dynare. If not, see . -global options_ oo_ M_ - -xparam = get_posterior_parameters(ptype); -gend = options_.nobs; +global options_ oo_ M_ bayestopt_ -% Read and demean data -rawdata = read_variables(options_.datafile,options_.varobs,[],options_.xls_sheet,options_.xls_range); -rawdata = rawdata(options_.first_obs:options_.first_obs+gend-1,:); -if options_.loglinear == 1 & ~options_.logdata - rawdata = log(rawdata); +if isfield(options_cond_fcst,'parameter_set') || isempty(options_cond_fcst.parameter_set) + options_cond_fcst.parameter_set = 'posterior_mode'; end -if options_.prefilter == 1 - bayestopt_.mean_varobs = mean(rawdata,1)'; - data = transpose(rawdata-repmat(bayestopt_.mean_varobs',gend,1)); + +if ischar(options_cond_fcst.parameter_set) + switch options_cond_fcst.parameter_set + case 'posterior_mode' + xparam = get_posterior_parameters('mode'); + case 'posterior_mean' + xparam = get_posterior_parameters('mean'); + case 'posterior_median' + xparam = get_posterior_parameters('median'); + case 'prior_mode' + xparam = bayestopt_.p5(:); + case 'prior_mean' + xparam = bayestopt_.p1; + otherwise + disp('eval_likelihood:: If the input argument is a string, then it has to be equal to:') + disp(' ''posterior_mode'', ') + disp(' ''posterior_mean'', ') + disp(' ''posterior_median'', ') + disp(' ''prior_mode'' or') + disp(' ''prior_mean''.') + error('imcforecast:: Wrong argument type!') + end else - data = transpose(rawdata); + xparam = options_cond_fcst.parameter_set; + if length(xparam)~=length(M_.params) + error('imcforecast:: The dimension of the vector of parameters doesn''t match the number of estimated parameters!') + end +end + +if isfield(options_cond_fcst,'replic') || isempty(options_cond_fcst.replic) + options_cond_fcst.replic = 5000; +end + +if isfield(options_cond_fcst,'periods') || isempty(options_cond_fcst.periods) + options_cond_fcst.replic = 40; +end + +if isfield(options_cond_fcst,'conf_sig') || isempty(options_cond_fcst.conf_sig) + options_cond_fcst.conf_sig = .8; end set_parameters(xparam); -[atT,innov,measurement_error,filtered_state_vector,ys,trend_coeff] = DsgeSmoother(xparam,gend,data,[],0); +n_varobs = size(options_.varobs,1); +rawdata = read_variables(options_.datafile,options_.varobs,[],options_.xls_sheet,options_.xls_range); +options_ = set_default_option(options_,'nobs',size(rawdata,1)-options_.first_obs+1); +gend = options_.nobs; +rawdata = rawdata(options_.first_obs:options_.first_obs+gend-1,:); +% Transform the data. +if options_.loglinear + if ~options_.logdata + rawdata = log(rawdata); + end +end +% Test if the data set is real. +if ~isreal(rawdata) + error('There are complex values in the data! Probably a wrong transformation') +end +% Detrend the data. +options_.missing_data = any(any(isnan(rawdata))); +if options_.prefilter == 1 + if options_.missing_data + bayestopt_.mean_varobs = zeros(n_varobs,1); + for variable=1:n_varobs + rdx = find(~isnan(rawdata(:,variable))); + m = mean(rawdata(rdx,variable)); + rawdata(rdx,variable) = rawdata(rdx,variable)-m; + bayestopt_.mean_varobs(variable) = m; + end + else + bayestopt_.mean_varobs = mean(rawdata,1)'; + rawdata = rawdata-repmat(bayestopt_.mean_varobs',gend,1); + end +end +data = transpose(rawdata); +% Handle the missing observations. +[data_index,number_of_observations,no_more_missing_observations] = describe_missing_data(data,gend,n_varobs); +missing_value = ~(number_of_observations == gend*n_varobs); -trend = repmat(ys,1,H+1); +[atT,innov,measurement_error,filtered_state_vector,ys,trend_coeff] = DsgeSmoother(xparam,gend,data,data_index,number_of_observations); + +trend = repmat(ys,1,options_cond_fcst.periods+1); for i=1:M_.endo_nbr j = strmatch(deblank(M_.endo_names(i,:)),options_.varobs,'exact'); if ~isempty(j) - trend(i,:) = trend(i,:)+trend_coeff(j)*(gend+(0:H)); + trend(i,:) = trend(i,:)+trend_coeff(j)*(gend+(0:options_cond_fcst.periods)); end end trend = trend(oo_.dr.order_var,:); @@ -54,30 +144,28 @@ InitState(:,1) = atT(:,end); sQ = sqrt(M_.Sigma_e); NumberOfStates = length(InitState); -FORCS1 = zeros(NumberOfStates,H+1,B); +FORCS1 = zeros(NumberOfStates,options_cond_fcst.periods+1,options_cond_fcst.replic); -for b=1:B +for b=1:options_cond_fcst.replic FORCS1(:,1,b) = InitState; end EndoSize = M_.endo_nbr; ExoSize = M_.exo_nbr; -n1 = size(cV,1); -n2 = size(cS,1); +n1 = size(constrained_vars,1); +n2 = size(options_cond_fcst.controlled_varexo,1); if n1 ~= n2 - disp('imcforecast :: Error!') - disp(['imcforecast :: The number of variables doesn''t match the number of shocks']) - return + error(['imcforecast:: The number of constrained variables doesn''t match the number of controlled shocks']) end idx = []; jdx = []; for i = 1:n1 - idx = [idx ; oo_.dr.inv_order_var(strmatch(deblank(cV(i,:)),M_.endo_names,'exact'))]; - jdx = [jdx ; strmatch(deblank(cS(i,:)),M_.exo_names,'exact')]; + idx = [idx ; oo_.dr.inv_order_var(strmatch(deblank(constrained_vars(i,:)),M_.endo_names,'exact'))]; + jdx = [jdx ; strmatch(deblank(options_cond_fcst.controlled_varexo(i,:)),M_.exo_names,'exact')]; end mv = zeros(n1,NumberOfStates); mu = zeros(ExoSize,n2); @@ -86,29 +174,28 @@ for i=1:n1 mu(jdx(i),i) = 1; end - -if (size(mcValue,2) == 1); - mcValue = mcValue*ones(1,cL); +if (size(constrained_paths,2) == 1); + constrained_paths = constrained_paths*ones(1,cL); else - cL = size(mcValue,2); + cL = size(constrained_paths,2); end randn('state',0); -for b=1:B - shocks = sQ*randn(ExoSize,H); - shocks(jdx,:) = zeros(length(jdx),H); - FORCS1(:,:,b) = mcforecast3(cL,H,mcValue,shocks,FORCS1(:,:,b),T,R,mv, mu)+trend; +for b=1:options_cond_fcst.replic + shocks = sQ*randn(ExoSize,options_cond_fcst.periods); + shocks(jdx,:) = zeros(length(jdx),options_cond_fcst.periods); + FORCS1(:,:,b) = mcforecast3(cL,options_cond_fcst.periods,constrained_paths,shocks,FORCS1(:,:,b),T,R,mv, mu)+trend; end mFORCS1 = mean(FORCS1,3); -tt = (1-ci)/2; -t1 = round(B*tt); -t2 = round(B*(1-tt)); +tt = (1-options_cond_fcst.conf_sig)/2; +t1 = round(options_cond_fcst.replic*tt); +t2 = round(options_cond_fcst.replic*(1-tt)); -forecasts.controled_variables = cV; -forecasts.instruments = cS; +forecasts.controled_variables = constrained_vars; +forecasts.instruments = options_cond_fcst.controlled_varexo; for i = 1:EndoSize eval(['forecasts.cond.mean.' deblank(M_.endo_names(oo_.dr.order_var(i),:)) ' = mFORCS1(i,:)'';']); @@ -119,22 +206,21 @@ end clear FORCS1; -FORCS2 = zeros(NumberOfStates,H+1,B); -for b=1:B +FORCS2 = zeros(NumberOfStates,options_cond_fcst.periods+1,options_cond_fcst.replic); +for b=1:options_cond_fcst.replic FORCS2(:,1,b) = InitState; end randn('state',0); -for b=1:B - shocks = sQ*randn(ExoSize,H); - shocks(jdx,:) = zeros(length(jdx),H); - FORCS2(:,:,b) = mcforecast3(0,H,mcValue,shocks,FORCS2(:,:,b),T,R,mv, mu)+trend; +for b=1:options_cond_fcst.replic + shocks = sQ*randn(ExoSize,options_cond_fcst.periods); + shocks(jdx,:) = zeros(length(jdx),options_cond_fcst.periods); + FORCS2(:,:,b) = mcforecast3(0,options_cond_fcst.periods,constrained_paths,shocks,FORCS2(:,:,b),T,R,mv, mu)+trend; end mFORCS2 = mean(FORCS2,3); - for i = 1:EndoSize eval(['forecasts.uncond.mean.' deblank(M_.endo_names(oo_.dr.order_var(i),:)) ' = mFORCS2(i,:)'';']); tmp = sort(squeeze(FORCS2(i,:,:))'); diff --git a/matlab/plot_icforecast.m b/matlab/plot_icforecast.m index 2b8fe5815..6dd76e4f4 100644 --- a/matlab/plot_icforecast.m +++ b/matlab/plot_icforecast.m @@ -1,6 +1,16 @@ -function plot_icforecast(Variable) +function plot_icforecast(Variables) +% Build plots for the conditional forecasts. +% +% INPUTS +% o Variables [char] m*x array holding the names of the endogenous variables to be plotted. +% +% OUTPUTS +% None. +% +% SPECIAL REQUIREMENTS +% This routine has to be called after imcforecast.m. -% Copyright (C) 2006 Dynare Team +% Copyright (C) 2006-2009 Dynare Team % % This file is part of Dynare. % @@ -19,30 +29,27 @@ function plot_icforecast(Variable) load conditional_forecasts; +for i=1:size(Variables,1) + eval(['ci1 = forecasts.cond.ci.' Variables(i,:) ';']) + eval(['m1 = forecasts.cond.mean.' Variables(i,:) ';']) + eval(['ci2 = forecasts.uncond.ci.' Variables(i,:) ';']) + eval(['m2 = forecasts.uncond.mean.' Variables(i,:) ';']) + build_figure(Variables(i,:),ci1,ci2,m1,m2); +end -eval(['ci1 = forecasts.cond.ci.' Variable ';']) -eval(['m1 = forecasts.cond.mean.' Variable ';']) -eval(['ci2 = forecasts.uncond.ci.' Variable ';']) -eval(['m2 = forecasts.uncond.mean.' Variable ';']) - - - -H = length(m1); - -% area(1:H,ci1(2,:),'FaceColor',[.9 .9 .9],'BaseValue',min([min(ci1(1,:)),min(ci2(1,:))])) - -h1 = area(1:H,ci1(2,1:H)) -set(h1,'BaseValue',min([min(ci1(1,:)),min(ci2(1,:))])) -set(h1,'FaceColor',[.9 .9 .9]) - -hold on -% area(1:H,ci1(1,:),'FaceColor',[1 1 1],'BaseValue',min([min(ci1(1,:)),min(ci2(1,:))])) -h2 = area(1:H,ci1(1,1:H)); -set(h2,'BaseValue',min([min(ci1(1,:)),min(ci2(1,:))])) -set(h2,'FaceColor',[1 1 1]) -plot(1:H,m1,'-k','linewidth',3) -plot(1:H,m2,'--k','linewidth',3) -plot(1:H,ci2(1,:),'--k','linewidth',1) -plot(1:H,ci2(2,:),'--k','linewidth',1) -axis tight -hold off \ No newline at end of file +function build_figure(name,ci1,ci2,m1,m2) + figure('Name',['Conditional forecast: ' name '.']); + H = length(m1); + h1 = area(1:H,ci1(2,1:H)); + set(h1,'BaseValue',min([min(ci1(1,:)),min(ci2(1,:))])) + set(h1,'FaceColor',[.9 .9 .9]) + hold on + h2 = area(1:H,ci1(1,1:H)); + set(h2,'BaseValue',min([min(ci1(1,:)),min(ci2(1,:))])) + set(h2,'FaceColor',[1 1 1]) + plot(1:H,m1,'-k','linewidth',3) + plot(1:H,m2,'--k','linewidth',3) + plot(1:H,ci2(1,:),'--k','linewidth',1) + plot(1:H,ci2(2,:),'--k','linewidth',1) + axis tight + hold off \ No newline at end of file