diff --git a/matlab/annualized_shock_decomposition.m b/matlab/annualized_shock_decomposition.m
new file mode 100644
index 000000000..e6f6479ef
--- /dev/null
+++ b/matlab/annualized_shock_decomposition.m
@@ -0,0 +1,195 @@
+function [z, endo_names, endo_names_tex, steady_state, i_var, oo_] = annualized_shock_decomposition(oo_, M_, opts, i_var, t0, t1, realtime_, vintage_, steady_state,GYTREND0,var_type,islog)
+% function oo_ = annualized_shock_decomposition(oo_,t0,options_.nobs);
+% Computes annualized shocks contribution to a simulated trajectory. The fields set are
+% oo_.annualized_shock_decomposition, oo_.annualized_realtime_shock_decomposition,
+% oo_.annualized_realtime_conditional_shock_decomposition and oo_.annualized_realtime_forecast_shock_decomposition.
+% Subfields are arrays n_var by nshock+2 by nperiods. The
+% first nshock columns store the respective shock contributions, column n+1
+% stores the role of the initial conditions, while column n+2 stores the
+% value of the smoothed variables. Both the variables and shocks are stored
+% in the order of endo_names and M_.exo_names, respectively.
+%
+% INPUTS
+% oo_: [structure] Storage of results
+% M_: [structure] Storage of model
+% opts: [structure] options for shock decomp
+% i_var: [array] index of vars
+% t0: [integer] first period
+% t1: [integer] last period
+% realtime_: [integer]
+% vintage_: [integer]
+% steady_state: [array] steady state value of quarterly (log-) level vars
+% GYTREND0: [array] growth of level of vars
+% var_type: [integer] flag for stock/flow/deflator
+% islog: [integer] flag for log-levels
+%
+% OUTPUTS
+% z: [matrix] shock decomp to plot
+% endo_names: [char] updated var names
+% endo_names_tex: [char] updated TeX var names
+% steady_state: [array] updated stady state of vars
+% i_var: [integer array] updated var indices to plot
+% oo_: [structure] Storage of results
+%
+% SPECIAL REQUIREMENTS
+% none
+
+% Copyright (C) 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 .
+
+nvar = length(i_var);
+
+% usual shock decomp
+z = oo_.shock_decomposition;
+z = z(i_var,:,:);
+if var_type==2,
+ gtxt = 'PHI'; % inflation rate
+ gtex = '\pi';
+else
+ gtxt = 'G'; % growth rate
+ gtex = 'g';
+end
+steady_state=steady_state(i_var);
+% endo_names = M_.endo_names(i_var,:);
+% endo_names_tex = M_.endo_names_tex(i_var,:);
+nterms = size(z,2);
+nfrcst = opts.forecast/4;
+
+for j=1:nvar
+ if j>1,
+ endo_names = char(endo_names,[deblank(M_.endo_names(i_var(j),:)) '_A']);
+ endo_names_tex = char(endo_names_tex,['{' deblank(M_.endo_names_tex(i_var(j),:)) '}^A']);
+ gendo_names = char(gendo_names,[gtxt endo_names(j,:)]);
+ gendo_names_tex = char(gendo_names_tex,[gtex '(' deblank(endo_names_tex(j,:)) ')']);
+ else
+ endo_names = [deblank(M_.endo_names(i_var(j),:)) '_A'];
+ endo_names_tex = ['{' deblank(M_.endo_names_tex(i_var(j),:)) '}^A'];
+ gendo_names = [gtxt endo_names(j,:)];
+ gendo_names_tex = [gtex '(' deblank(endo_names_tex(j,:)) ')'];
+ end
+ for k =1:nterms,
+ [za(j,k,:), steady_state_a(j,1), gza(j,k,:), steady_state_ga(j,1)] = ...
+ quarterly2annual(squeeze(z(j,k,t0:end)),steady_state(j),GYTREND0,var_type,islog);
+ end
+ ztmp=squeeze(za(j,:,:));
+ zscale = sum(ztmp(1:end-1,:))./ztmp(end,:);
+ ztmp(1:end-1,:) = ztmp(1:end-1,:)./repmat(zscale,[nterms-1,1]);
+ gztmp=squeeze(gza(j,:,:));
+ gscale = sum(gztmp(1:end-1,:))./ gztmp(end,:);
+ gztmp(1:end-1,:) = gztmp(1:end-1,:)./repmat(gscale,[nterms-1,1]);
+ za(j,:,:) = ztmp;
+ gza(j,:,:) = gztmp;
+end
+
+z=cat(1,za,gza);
+oo_.annualized_shock_decomposition=z;
+endo_names = char(endo_names,gendo_names);
+endo_names_tex = char(endo_names_tex,gendo_names_tex);
+
+% realtime
+init=1;
+for i=t0:4:t1,
+ yr=floor(i/4);
+ za=[];
+ gza=[];
+ z = oo_.realtime_shock_decomposition.(['time_' int2str(i)]);
+ z = z(i_var,:,:);
+
+ for j=1:nvar
+ for k =nterms:-1:1,
+% if knfrcst
+ oo_.annualized_realtime_conditional_shock_decomposition.(['yr_' int2str(yr-nfrcst)]) = ...
+ oo_.annualized_realtime_shock_decomposition.pool(:,:,yr-nfrcst:end) - ...
+ oo_.annualized_realtime_forecast_shock_decomposition.(['yr_' int2str(yr-nfrcst)]);
+ oo_.annualized_realtime_conditional_shock_decomposition.(['yr_' int2str(yr-nfrcst)])(:,end-1,:) = ...
+ oo_.annualized_realtime_forecast_shock_decomposition.(['yr_' int2str(yr-nfrcst)])(:,end,:);
+ end
+ end
+% ztmp=oo_.realtime_shock_decomposition.pool(:,:,21:29)-oo_.realtime_forecast_shock_decomposition.time_21;
+
+
+
+ init=init+1;
+end
+
+i_var=1:2*nvar;
+steady_state = [steady_state_a;steady_state_ga];
+
+
+switch realtime_
+
+ case 0
+ z = oo_.annualized_shock_decomposition;
+
+ case 1 % realtime
+ if vintage_
+ z = oo_.annualized_realtime_shock_decomposition.(['yr_' int2str(floor(vintage_/4))]);
+ else
+ z = oo_.annualized_realtime_shock_decomposition.pool;
+ end
+
+ case 2 % conditional
+ if vintage_
+ z = oo_.annualized_realtime_conditional_shock_decomposition.(['yr_' int2str(floor(vintage_/4))]);
+ else
+ error();
+ end
+
+ case 3 % forecast
+ if vintage_
+ z = oo_.annualized_realtime_forecast_shock_decomposition.(['yr_' int2str(floor(vintage_/4))]);
+ else
+ error()
+ end
+end
diff --git a/matlab/plot_shock_decomposition.m b/matlab/plot_shock_decomposition.m
index 6cdd3a23d..5935aa3f6 100644
--- a/matlab/plot_shock_decomposition.m
+++ b/matlab/plot_shock_decomposition.m
@@ -67,6 +67,12 @@ end
write_xls = options_.shock_decomp.write_xls;
initial_date = options_.initial_date;
+
+if isfield(options_.shock_decomp,'q2a'), % private trap for aoa calls
+ q2a=options_.shock_decomp.q2a;
+else
+ q2a=0;
+end
switch realtime_
@@ -85,7 +91,7 @@ switch realtime_
case 2 % conditional
if vintage_
- z = oo_.conditional_shock_decomposition.(['time_' int2str(vintage_)]);
+ z = oo_.realtime_conditional_shock_decomposition.(['time_' int2str(vintage_)]);
initial_date = options_.initial_date+vintage_-1;
fig_names1=[fig_names ' ' int2str(forecast_) '-step ahead conditional forecast (given ' char(initial_date) ')'];
else
@@ -170,14 +176,47 @@ switch type
case 'aoa'
if isempty(initial_date),
- initial_date = dates('1Y');
t0=4;
+ initial_date = dates('1Y');
+ else
+ initial_date0 = dates([int2str(initial_date.time(1)) 'Y']);
+ if initial_date.time(2)==1,
+ t0=4;
+ else
+ t0=4+(4-initial_date.time(2)+1);
+ initial_date1=initial_date0+1;
+ end
+ end
+ if q2a
+ var_type=1;
+ islog=0;
+ GYTREND0 = 0;
+ if isfield(options_.shock_decomp,'var_type'), % private trap for aoa calls
+ var_type=options_.shock_decomp.var_type;
+ end
+ if isfield(options_.shock_decomp,'islog'), % private trap for aoa calls
+ islog=options_.shock_decomp.islog;
+ end
+ if isfield(options_.shock_decomp,'GYTREND0'), % private trap for aoa calls
+ GYTREND0=options_.shock_decomp.GYTREND0;
+ end
+
+
+ [za, endo_names, endo_names_tex, steady_state, i_var, oo_] = ...
+ annualized_shock_decomposition(oo_,M_, options_.shock_decomp, i_var, t0, options_.nobs, realtime_, vintage_, steady_state,GYTREND0,var_type,islog);
+ if realtime_<2
+ initial_date = initial_date1;
+ else
+ initial_date = initial_date0;
+ end
+ z = za;
+ M_.endo_names = endo_names;
+ M_.endo_names_tex = endo_names_tex;
else
t0=4-initial_date.time(2)+1;
- initial_date = dates([int2str(initial_date.time(1)) 'Y']);
+ initial_date = initial_date0;
+ z=z(:,:,t0:4:end);
end
- z=z(:,:,t0:4:end);
-
otherwise
error('plot_shock_decomposition:: Wrong type')
diff --git a/matlab/utilities/dataset/lagged.m b/matlab/utilities/dataset/lagged.m
new file mode 100644
index 000000000..31a58909a
--- /dev/null
+++ b/matlab/utilities/dataset/lagged.m
@@ -0,0 +1,32 @@
+function xlag = lagged(x, n);
+% xlag = lagged(x, n);
+% applies n-lags backward shift operator to x
+%
+% INPUTS
+% x = time series
+% n = number of backward shifts [DEFAULT=1]
+%
+% OUTPUT
+% xlag = backward shifted series
+
+% Copyright (C) 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 .
+
+if nargin==1, n=1; end
+
+x=x(:);
+xlag=[NaN(n,1); x(1:end-n)];
diff --git a/matlab/utilities/dataset/quarterly2annual.m b/matlab/utilities/dataset/quarterly2annual.m
new file mode 100644
index 000000000..d615e923f
--- /dev/null
+++ b/matlab/utilities/dataset/quarterly2annual.m
@@ -0,0 +1,79 @@
+function [ya, yass, gya, gyass] = quarterly2annual(y,yss,GYTREND0,type,islog)
+% function [ya, yass, gya, gyass] = quarterly2annual(y,yss,GYTREND0,type,islog)
+% transforms quarterly (log-)level time series to annual level and growth rate
+% it accounts for stock/flow/deflator series.
+%
+% INPUTS
+% y quarterly time series
+% yss steady state of y
+% GYTREND0 growth rate of y
+% type 1 flow (default)
+% 2 deflator
+% 0 stock
+% islog 0 level (default)
+% 1 log-level
+%
+% OUTPUTS
+% ya annual (log-)level
+% yass annual steadystate (log-)level
+% gya annual growth rate
+% gyass annual growth rate steadystate
+
+% Copyright (C) 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 .
+
+if nargin ==0,
+ disp('[ya, yass, gya, gyass] = quarterly2annual(y,yss,GYTREND0,type,islog);')
+ return
+end
+
+if nargin<4 || isempty(type),
+ type=1;
+end
+if nargin<5 || isempty(islog),
+ islog=0;
+end
+if islog
+ y=exp(y+yss);
+ yss=exp(yss);
+ y=y-yss;
+end
+switch type
+ case 1
+ yass = yss*(exp(-GYTREND0*3)+exp(-GYTREND0*2)+exp(-GYTREND0)+1);
+ tmp = lagged(y,3)*exp(-GYTREND0*3)+lagged(y,2)*exp(-GYTREND0*2)+lagged(y,1)*exp(-GYTREND0)+y; % annualized level
+ case 2
+ yass = yss*(exp(-GYTREND0*3)+exp(-GYTREND0*2)+exp(-GYTREND0)+1)/4;
+ tmp = (lagged(y,3)*exp(-GYTREND0*3)+lagged(y,2)*exp(-GYTREND0*2)+lagged(y,1)*exp(-GYTREND0)+y)/4; % annualized level
+ case 3
+ yass=yss;
+ tmp = y;
+ otherwise
+ error('Wrong type input')
+end
+
+ya = tmp(4:4:end);
+% annual growth rate
+gyass = GYTREND0*4;
+gya = (ya+yass)./(lagged(ya,1)+yass).*exp(4*GYTREND0)-1-gyass;
+
+if islog
+ ya=log(ya+yass);
+ yass=log(yass);
+ ya=ya-yass;
+end
+