diff --git a/matlab/+get/irf.m b/matlab/+get/irf.m
new file mode 100755
index 000000000..e247d2054
--- /dev/null
+++ b/matlab/+get/irf.m
@@ -0,0 +1,47 @@
+
+function y0 = get_irf(exo,varargin)
+% function x = get_irf(exoname, vname1, vname2, ...)
+% returns IRF to individual exogenous for a list of variables and adds the
+% steady state
+%
+% INPUTS:
+%   exo: exo variable name
+%   vname1, vname2, ... :  list of variable names
+%
+% OUTPUTS
+%   x:      irf matrix [time x number of variables]
+%
+% SPECIAL REQUIREMENTS
+%   none
+
+% Copyright (C) 2019 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 <http://www.gnu.org/licenses/>.
+
+global M_ oo_
+
+ys_ = [oo_.steady_state];
+y0=zeros(length(oo_.irfs.([varargin{1} '_' exo]))+1,length(varargin));
+
+
+[i_var,nvar] = varlist_indices(varargin,M_.endo_names);
+
+
+for j=1:nvar
+%     mfys = strmatch(varargin{j},lgy_,'exact');
+    y0(:,j)=[0; oo_.irfs.([ varargin{j} '_' exo ])']+ys_(i_var(j));
+end
+
diff --git a/matlab/+get/mean.m b/matlab/+get/mean.m
new file mode 100755
index 000000000..f307210dc
--- /dev/null
+++ b/matlab/+get/mean.m
@@ -0,0 +1,56 @@
+function y0 = get_mean(varargin)
+% function x = get_mean(vname1, vname2, <order>)
+% returns the steady-state of a variable identified by its name
+%
+% INPUTS:
+%   vname1, vname2, ... :  list of variable names
+%   order: if integer 1 or 2, optionally last input can trigger the order
+%   at which steady state is computed
+%
+% OUTPUTS
+%   x:      steady state values
+%
+% SPECIAL REQUIREMENTS
+%   none
+
+% Copyright (C) 2019 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 <http://www.gnu.org/licenses/>.
+
+global M_ oo_ options_ 
+
+if ~isempty(regexp(varargin{end},'\d','ONCE')) && isempty(regexp(varargin{end},'\D','ONCE'))
+    order=eval(varargin{end});
+else
+    order=1;
+end
+if order==1
+     ys_ = oo_.steady_state;
+     ys_ = evaluate_steady_state(ys_,M_,options_,oo_,1);
+elseif order==2
+    ys_ = oo_.dr.ys;
+    ys_(oo_.dr.order_var)=ys_(oo_.dr.order_var)+oo_.dr.ghs2./2;
+else
+   return
+end
+lgy_ = M_.endo_names;
+
+mfys=nan(length(varargin),1);
+for j=1:length(varargin)
+    mfys(j) = find(strcmp(varargin{j},lgy_));
+end
+
+y0 = ys_(mfys);
diff --git a/matlab/+get/shock_stderr_by_name.m b/matlab/+get/shock_stderr_by_name.m
new file mode 100755
index 000000000..b1143a69a
--- /dev/null
+++ b/matlab/+get/shock_stderr_by_name.m
@@ -0,0 +1,39 @@
+function x = get_shock_stderr_by_name(exoname)
+% function x = get_shock_stderr_by_name(exoname)
+% returns the value of a shock identified by its name
+%  
+% INPUTS:
+%   exoname:  shock name
+%
+% OUTPUTS
+%   x:      shock value
+%
+% SPECIAL REQUIREMENTS
+%   none
+
+% Copyright (C) 2019 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 <http://www.gnu.org/licenses/>.
+
+global M_
+
+i = find(strcmp(exoname,M_.exo_names));
+
+if isempty(i)
+    error(['Can''t find shock ', exoname])
+end
+
+x = sqrt(M_.Sigma_e(i,i));
diff --git a/matlab/+get/smooth.m b/matlab/+get/smooth.m
new file mode 100755
index 000000000..c81cbbc1b
--- /dev/null
+++ b/matlab/+get/smooth.m
@@ -0,0 +1,46 @@
+
+function y0 = get_smooth(varargin)
+% function x = get_smooth(vname1, vname2, )
+% returns smoothed variables or shocks identified by their name
+%
+% INPUTS:
+%   vname1, vname2, ... :  list of variable/shock names
+%
+% OUTPUTS
+%   x:      smoothed variables [T x number of variables]
+%
+% SPECIAL REQUIREMENTS
+%   none
+
+% Copyright (C) 2019 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 <http://www.gnu.org/licenses/>.
+global oo_
+
+SmoothedVariables=[struct2cell(oo_.SmoothedVariables); struct2cell(oo_.SmoothedShocks)];
+my_field_names = [fieldnames(oo_.SmoothedVariables); fieldnames(oo_.SmoothedShocks)];
+isvar=zeros(length(SmoothedVariables),1);
+for jf = 1:length(SmoothedVariables)
+    isvar(jf)=~(isstruct(SmoothedVariables{jf}));
+end
+SmoothedVariables=cell2struct(SmoothedVariables(logical(isvar)),my_field_names(logical(isvar)));
+
+
+y0=zeros(length(SmoothedVariables.(varargin{1})),length(varargin));
+for j=1:length(varargin)
+    y0(:,j)=SmoothedVariables.(varargin{j}); 
+end
+
diff --git a/matlab/+get/update.m b/matlab/+get/update.m
new file mode 100755
index 000000000..994a6e47b
--- /dev/null
+++ b/matlab/+get/update.m
@@ -0,0 +1,36 @@
+function y0 = get_update(varargin)
+% function x = get_update(vname1, vname2, )
+% returns updated variables identified by their name
+%
+% INPUTS:
+%   vname1, vname2, ... :  list of variable names
+%
+% OUTPUTS
+%   x:      smoothed variables [T x number of variables]
+%
+% SPECIAL REQUIREMENTS
+%   none
+
+% Copyright (C) 2019 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 <http://www.gnu.org/licenses/>.
+global oo_ 
+
+y0=zeros(length(oo_.UpdatedVariables.(varargin{1})),length(varargin));
+for j=1:length(varargin)
+    y0(:,j)=oo_.UpdatedVariables.(varargin{j}); 
+end
+
diff --git a/matlab/WriteShockDecomp2Excel.m b/matlab/WriteShockDecomp2Excel.m
index f9c18fc7a..a11110455 100644
--- a/matlab/WriteShockDecomp2Excel.m
+++ b/matlab/WriteShockDecomp2Excel.m
@@ -100,8 +100,8 @@ for j=1:nvar
         comp_nbr=18;
     end
 
-    d0(1,:)=[{'Decomposition'} cellstr(labels(1:comp_nbr,:))' {'Smoot Var'}];
-    d0=[d0; num2cell([x' z1'])];
+    d0(1,:)=[{'Decomposition'} cellstr(labels(1:comp_nbr,:))' {'Smoot Var'} {'Steady State'}];
+    d0=[d0; num2cell([x' z1' ]), [num2cell(SteadyState(i_var(j))); cell(size(z1,2)-1,1)]];
     LastRow=size(d0,1);
     if use_shock_groups
         d0(LastRow+2,1)={'Legend.'};
diff --git a/matlab/annualized_shock_decomposition.m b/matlab/annualized_shock_decomposition.m
index e7a645764..03b5feb44 100644
--- a/matlab/annualized_shock_decomposition.m
+++ b/matlab/annualized_shock_decomposition.m
@@ -158,7 +158,7 @@ if realtime_ && isstruct(oo_) && isfield(oo_, 'realtime_shock_decomposition')
         myopts.plot_shock_decomp.realtime=1;
         myopts.plot_shock_decomp.vintage=i;
         % retrieve quarterly shock decomp
-        z = plot_shock_decomposition(M_,oo_,myopts,[]);
+        [z, ~] = plot_shock_decomposition(M_,oo_,myopts,[]);
         zdim = size(z);
         z = z(i_var,:,:);
         if isstruct(aux)
@@ -185,13 +185,14 @@ if realtime_ && isstruct(oo_) && isfield(oo_, 'realtime_shock_decomposition')
             if qvintage_>i-4 && qvintage_<i
                 myopts.plot_shock_decomp.vintage=qvintage_;
                 % retrieve quarterly shock decomp
-                z = plot_shock_decomposition(M_,oo_,myopts,[]);
+                [z, ~] = plot_shock_decomposition(M_,oo_,myopts,[]);
                 z(:,:,end+1:zdim(3))=nan; % fill with nan's remaining time points to reach Q4
                 z = z(i_var,:,:);
                 if isstruct(aux)
                     if ischar(aux0.y)
                         % retrieve quarterly shock decomp for aux variable
                         [y_aux, steady_state_aux] = plot_shock_decomposition(M_,oo_,myopts,aux0.y);
+                        y_aux(:,:,end+1:zdim(3))=nan; % fill with nan's remaining time points to reach Q4
                         aux.y=y_aux;
                         aux.yss=steady_state_aux;
                     end
@@ -202,6 +203,7 @@ if realtime_ && isstruct(oo_) && isfield(oo_, 'realtime_shock_decomposition')
 
             end
             oo_.annualized_realtime_forecast_shock_decomposition.(['yr_' int2str(yr)]) = z(:,:,end-nfrcst:end);
+            oo_.annualized_realtime_forecast_shock_decomposition.pool(:,:,yr+1) = squeeze(z(:,:,end-nfrcst+1));
             if init>nfrcst
                 oo_.annualized_realtime_conditional_shock_decomposition.(['yr_' int2str(yr-nfrcst)]) = ...
                     oo_.annualized_realtime_shock_decomposition.pool(:,:,yr-nfrcst:end) - ...
@@ -223,8 +225,12 @@ if realtime_ && isstruct(oo_) && isfield(oo_, 'realtime_shock_decomposition')
                             oo_.annualized_realtime_forecast_shock_decomposition.(['yr_' int2str(yr-my_forecast_)])(:,end,1:my_forecast_+1);
                         oo_.annualized_realtime_conditional_shock_decomposition.(['yr_' int2str(yr-my_forecast_)])(:,end,:) = ...
                             oo_.annualized_realtime_shock_decomposition.pool(:,end,yr-my_forecast_:yr);
+                        oo_.annualized_realtime_conditional_shock_decomposition.pool(:,:,yr-my_forecast_+1) = ...
+                            oo_.annualized_realtime_conditional_shock_decomposition.(['yr_' int2str(yr-my_forecast_)])(:,:,2);
                     end
                 end
+                oo_.annualized_realtime_conditional_shock_decomposition.pool(:,:,yr-nfrcst+1) = ...
+                    oo_.annualized_realtime_conditional_shock_decomposition.(['yr_' int2str(yr-nfrcst)])(:,:,2);
             end
         end
         % ztmp=oo_.realtime_shock_decomposition.pool(:,:,21:29)-oo_.realtime_forecast_shock_decomposition.time_21;
@@ -251,14 +257,14 @@ if realtime_ && isstruct(oo_) && isfield(oo_, 'realtime_shock_decomposition')
         if vintage_
             z = oo_.annualized_realtime_conditional_shock_decomposition.(['yr_' int2str(floor(vintage_/4))]);
         else
-            error();
+            z = oo_.annualized_realtime_conditional_shock_decomposition.pool;
         end
 
       case 3 % forecast
         if vintage_
             z = oo_.annualized_realtime_forecast_shock_decomposition.(['yr_' int2str(floor(vintage_/4))]);
         else
-            error()
+            z = oo_.annualized_realtime_forecast_shock_decomposition.pool;
         end
     end
 end
diff --git a/matlab/expand_group.m b/matlab/expand_group.m
index 8011b471b..0992e4967 100644
--- a/matlab/expand_group.m
+++ b/matlab/expand_group.m
@@ -67,4 +67,4 @@ options.plot_shock_decomp.colormap = MAP;
 options.plot_shock_decomp.endo_names = mydata.endo_names;
 options.plot_shock_decomp.endo_names_tex = mydata.endo_names_tex;
 
-plot_shock_decomposition(M,oo,options,var_list_);
+plot_shock_decomposition(M,oo,options,{var_list_});
diff --git a/matlab/graph_decomp.m b/matlab/graph_decomp.m
index b0c865f13..b3190bf7f 100644
--- a/matlab/graph_decomp.m
+++ b/matlab/graph_decomp.m
@@ -148,6 +148,11 @@ for j=1:nvar
     ax=axes('Position',[0.1 0.1 0.6 0.8],'box','on');
     %     plot(ax,x(2:end),z1(end,:),'k-','LineWidth',2)
     %     axis(ax,[xmin xmax ymin ymax]);
+    if strcmp('aoa',DynareOptions.plot_shock_decomp.type)
+        bgap = 0.15;
+    else
+        bgap = 0;
+    end        
     hold on;
     for i=1:gend
         i_1 = i-1;
@@ -156,10 +161,10 @@ for j=1:nvar
         for k = 1:comp_nbr
             zz = z1(k,i);
             if zz > 0
-                fill([x(i) x(i) x(i+1) x(i+1)]+(1/freq/2),[yp yp+zz yp+zz yp],k);
+                fill([x(i)+bgap x(i)+bgap x(i+1)-bgap x(i+1)-bgap]+(1/freq/2),[yp yp+zz yp+zz yp],k);
                 yp = yp+zz;
             else
-                fill([x(i) x(i) x(i+1) x(i+1)]+(1/freq/2),[ym ym+zz ym+zz ym],k);
+                fill([x(i)+bgap x(i)+bgap x(i+1)-bgap x(i+1)-bgap]+(1/freq/2),[ym ym+zz ym+zz ym],k);
                 ym = ym+zz;
             end
             hold on;
@@ -220,7 +225,7 @@ for j=1:nvar
                 c = uicontextmenu;
                 hl.UIContextMenu=c;
                 browse_menu = uimenu(c,'Label','Browse group');
-                expand_menu = uimenu(c,'Label','Expand group','Callback',['expand_group(''' mydata.plot_shock_decomp.use_shock_groups ''',''' deblank(mydata.plot_shock_decomp.orig_varlist{j}) ''',' int2str(i) ')']);
+                expand_menu = uimenu(c,'Label','Expand group','Callback',['expand_group(''' mydata.plot_shock_decomp.use_shock_groups ''',''' mydata.plot_shock_decomp.orig_varlist{j} ''',' int2str(i) ')']);
                 set(expand_menu,'UserData',mydata,'Tag',['group' int2str(i)]);
                 for jmember = mydata.shock_group.shocks
                     uimenu('parent',browse_menu,'Label',char(jmember))
@@ -252,6 +257,7 @@ for j=1:nvar
             dyn_saveas(fhandle,[DynareOptions.plot_shock_decomp.filepath, filesep, DynareModel.fname,preamble_figname,endo_names{i_var(j)},fig_mode1,fig_name],DynareOptions.plot_shock_decomp.nodisplay,DynareOptions.plot_shock_decomp.graph_format);
         end
     end
+    
 end
 
 %% write LaTeX-Footer
diff --git a/matlab/graph_decomp_detail.m b/matlab/graph_decomp_detail.m
index 7ff2be01b..448abae6c 100644
--- a/matlab/graph_decomp_detail.m
+++ b/matlab/graph_decomp_detail.m
@@ -24,7 +24,7 @@ function []=graph_decomp_detail(z,shock_names,endo_names,i_var,initial_date,Dyna
 % 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.
-%
+%xf
 % You should have received a copy of the GNU General Public License
 % along with Dynare.  If not, see <http://www.gnu.org/licenses/>.
 
@@ -36,6 +36,7 @@ fig_mode='';
 fig_mode1='';
 % fig_name='';
 % screen_shocks=0;
+initval = DynareOptions.plot_shock_decomp.initval;
 use_shock_groups = DynareOptions.plot_shock_decomp.use_shock_groups;
 if use_shock_groups
     shock_groups = DynareModel.shock_groups.(use_shock_groups);
@@ -59,6 +60,11 @@ if isfield(opts_decomp,'init_cond_decomp')
 else
     init_cond_decomp = 0;
 end
+if isfield(opts_decomp,'min_nrows')
+    min_nrows = opts_decomp.min_nrows ;
+else
+    max_nrows = 6;
+end
 screen_shocks = opts_decomp.screen_shocks;
 if ~isempty(DynareOptions.plot_shock_decomp.use_shock_groups) || comp_nbr<=18
     screen_shocks=0;
@@ -138,9 +144,13 @@ end
 ncol=3;
 nrow=ceil(comp_nbr/ncol);
 ntotrow = nrow;
-nrow = min(ntotrow, 6);
+nrow = min(ntotrow, max_nrows);
 nfigs = ceil(ntotrow/nrow);
-labels = char(char(shock_names),'Initial values');
+if initval
+    labels = char(char(shock_names),'All shocks');
+else
+    labels = char(char(shock_names),'Initial values');
+end
 if ~(screen_shocks && comp_nbr>18)
     screen_shocks=0;
 end
@@ -150,7 +160,11 @@ for j=1:nvar
     z1 = squeeze(z(i_var(j),:,:));
     if screen_shocks,
         [~, isort] = sort(mean(abs(z1(1:end-2,:)')), 'descend');
-        labels = char(char(shock_names(isort(1:16))),'Others', 'Initial values');
+        if initval
+            labels = char(char(shock_names(isort(1:16))),'Others', 'All shocks');
+        else
+            labels = char(char(shock_names(isort(1:16))),'Others', 'Initial values');
+        end
         zres = sum(z1(isort(17:end),:),1);
         z1 = [z1(isort(1:16),:); zres; z1(comp_nbr0:end,:)];
         comp_nbr=18;
@@ -221,7 +235,7 @@ for j=1:nvar
                     c = uicontextmenu;
                     hax.UIContextMenu=c;
                     browse_menu = uimenu(c,'Label','Browse group');
-                    expand_menu = uimenu(c,'Label','Expand group','Callback',['expand_group(''' mydata.plot_shock_decomp.use_shock_groups ''',''' deblank(mydata.plot_shock_decomp.orig_varlist{j}) ''',' int2str(ic) ')']);
+                    expand_menu = uimenu(c,'Label','Expand group','Callback',['expand_group(''' mydata.plot_shock_decomp.use_shock_groups ''',''' mydata.plot_shock_decomp.orig_varlist{j} ''',' int2str(ic) ')']);
                     set(expand_menu,'UserData',mydata,'Tag',['group' int2str(ic)]);
                     for jmember = mydata.shock_group.shocks
                         uimenu('parent',browse_menu,'Label',char(jmember))
diff --git a/matlab/initial_condition_decomposition.m b/matlab/initial_condition_decomposition.m
index f6164cab2..12e39691a 100644
--- a/matlab/initial_condition_decomposition.m
+++ b/matlab/initial_condition_decomposition.m
@@ -62,8 +62,10 @@ if isempty(varlist)
     varlist = M_.endo_names(1:M_.orig_endo_nbr);
 end
 
-[i_var, nvar, index_uniques] = varlist_indices(varlist, M_.endo_names);
-varlist = varlist(index_uniques);
+if ~isequal(varlist,0)
+    [i_var, nvar, index_uniques] = varlist_indices(varlist, M_.endo_names);
+    varlist = varlist(index_uniques);
+end
 
 % number of variables
 endo_nbr = M_.endo_nbr;
@@ -83,7 +85,16 @@ if isempty(parameter_set)
     end
 end
 
-if ~isfield(oo_,'initval_decomposition')
+if ~isfield(oo_,'initval_decomposition') || isequal(varlist,0)
+    if isfield(oo_,'shock_decomposition_info') && isfield(oo_.shock_decomposition_info,'i_var')
+        if isfield (oo_,'realtime_conditional_shock_decomposition') ...
+                || isfield (oo_,'realtime_forecast_shock_decomposition') ...
+                || isfield (oo_,'realtime_shock_decomposition') ...
+                || isfield (oo_,'conditional_shock_decomposition') ...
+                || isfield (oo_,'shock_decomposition')
+            error('initval_decomposition::squeezed shock decompositions are already stored in oo_')
+        end
+    end
     options_.selected_variables_only = 0; %make sure all variables are stored
     options_.plot_priors=0;
     [oo,M,~,~,Smoothed_Variables_deviation_from_mean] = evaluate_smoother(parameter_set,varlist,M_,oo_,options_,bayestopt_,estim_params_);
@@ -129,27 +140,30 @@ if ~isfield(oo_,'initval_decomposition')
 
     end
 
-
     oo_.initval_decomposition = z;
 end
-% if ~options_.no_graph.shock_decomposition
-oo=oo_;
-oo.shock_decomposition = oo_.initval_decomposition;
-if ~isempty(init2shocks)
-    init2shocks = M_.init2shocks.(init2shocks);
-    n=size(init2shocks,1);
-    for i=1:n
-        j=strmatch(init2shocks{i}{1},M_.endo_names,'exact');
-        oo.shock_decomposition(:,end-1,:)=oo.shock_decomposition(:,j,:)+oo.shock_decomposition(:,end-1,:);
-        oo.shock_decomposition(:,j,:)=0;
-    end    
-end    
-M_.exo_names = M_.endo_names;
-M_.exo_nbr = M_.endo_nbr;
-options_.plot_shock_decomp.realtime=0;
-options_.plot_shock_decomp.screen_shocks=1;
-options_.plot_shock_decomp.use_shock_groups = '';
-options_.plot_shock_decomp.init_cond_decomp = 1; % private flag to plotting utilities
 
-plot_shock_decomposition(M_,oo,options_,varlist);
-% end
+% when varlist==0, we only store results in oo_ and do not make any plot
+if ~isequal(varlist,0)
+    
+    % if ~options_.no_graph.shock_decomposition
+    oo=oo_;
+    oo.shock_decomposition = oo_.initval_decomposition;
+    if ~isempty(init2shocks)
+        init2shocks = M_.init2shocks.(init2shocks);
+        n=size(init2shocks,1);
+        for i=1:n
+            j=strmatch(init2shocks{i}{1},M_.endo_names,'exact');
+            oo.shock_decomposition(:,end-1,:)=oo.shock_decomposition(:,j,:)+oo.shock_decomposition(:,end-1,:);
+            oo.shock_decomposition(:,j,:)=0;
+        end
+    end
+    M_.exo_names = M_.endo_names;
+    M_.exo_nbr = M_.endo_nbr;
+    options_.plot_shock_decomp.realtime=0;
+    options_.plot_shock_decomp.screen_shocks=1;
+    options_.plot_shock_decomp.use_shock_groups = '';
+    options_.plot_shock_decomp.init_cond_decomp = 1; % private flag to plotting utilities
+    
+    plot_shock_decomposition(M_,oo,options_,varlist);
+end
diff --git a/matlab/plot_shock_decomposition.m b/matlab/plot_shock_decomposition.m
index 3033b2e22..50b978cbc 100644
--- a/matlab/plot_shock_decomposition.m
+++ b/matlab/plot_shock_decomposition.m
@@ -1,4 +1,4 @@
-function [z, steady_state] = plot_shock_decomposition(M_,oo_,options_,varlist)
+function [out, steady_state] = plot_shock_decomposition(M_,oo_,options_,varlist)
 % function plot_shock_decomposition(M_,oo_,options_,varlist)
 % Plots the results of shock_decomposition
 %
@@ -31,9 +31,23 @@ function [z, steady_state] = plot_shock_decomposition(M_,oo_,options_,varlist)
 options_.nodisplay = options_.plot_shock_decomp.nodisplay;
 options_.graph_format = options_.plot_shock_decomp.graph_format;
 
+if ~isfield(oo_,'shock_decomposition_info')
+    oo_.shock_decomposition_info = struct();
+end
+if ~isfield(oo_,'plot_shock_decomposition_info')
+    oo_.plot_shock_decomposition_info = struct();
+end
+
+out=oo_;
 % indices of endogenous variables
-if isempty(varlist)
-    varlist = M_.endo_names(1:M_.orig_endo_nbr);
+exist_varlist = 1;
+if size(varlist,1) == 0
+    exist_varlist = 0;
+    if size( M_.endo_names,1) >= M_.orig_endo_nbr
+        varlist = M_.endo_names(1:M_.orig_endo_nbr);
+    else
+        varlist = M_.endo_names;
+    end
 end
 
 if isfield(options_.plot_shock_decomp,'init2shocks') % private trap for uimenu calls
@@ -45,10 +59,66 @@ if ~isempty(init2shocks)
     init2shocks=M_.init2shocks.(init2shocks);
 end
 
+if isfield(oo_.shock_decomposition_info,'i_var') && (M_.endo_nbr>=M_.orig_endo_nbr)
+    M_.endo_names = M_.endo_names(oo_.shock_decomposition_info.i_var,:);
+    M_.endo_names_tex = M_.endo_names_tex(oo_.shock_decomposition_info.i_var,:);
+    M_.endo_nbr = length( oo_.shock_decomposition_info.i_var );
+end
+try
+    [i_var,nvar,index_uniques] = varlist_indices(varlist,M_.endo_names);
+catch ME
+    if isfield(oo_.shock_decomposition_info,'i_var')
+        warning('shock decomp results for some input variable was not stored: I recompute all decompositions')
+        M_ = evalin('base','M_');
+        bayestopt_ = evalin('base','bayestopt_');
+        estim_params_ = evalin('base','estim_params_');
+        options_.no_graph.shock_decomposition=1; % force nograph in computing decompositions!
+        oo_.shock_decomposition_info = rmfield(oo_.shock_decomposition_info,'i_var');
+        var_list_ = char();
+        disp('recomputing shock decomposition ...')
+        [oo_,M_]= shock_decomposition(M_,oo_,options_,var_list_,bayestopt_,estim_params_);
+        if isfield(oo_,'realtime_shock_decomposition') || options_.plot_shock_decomp.realtime
+            disp('recomputing realtime shock decomposition ...')
+            oo_ = realtime_shock_decomposition(M_,oo_,options_,var_list_,bayestopt_,estim_params_);
+        end
+        if isfield(oo_,'initval_decomposition')
+            disp('recomputing initval shock decomposition ...')
+            oo_ = initial_condition_decomposition(M_,oo_,options_,0,bayestopt_,estim_params_);   
+        end
+        [i_var,nvar,index_uniques] = varlist_indices(varlist,M_.endo_names);
+        out = oo_;
+    else
+        rethrow(ME)
+    end
+end
 
-[i_var, ~, index_uniques] = varlist_indices(varlist, M_.endo_names);
 varlist = varlist(index_uniques);
 
+if ~isfield(out.shock_decomposition_info,'i_var') && exist_varlist
+    if ~isfield(out.plot_shock_decomposition_info,'i_var')
+        out.plot_shock_decomposition_info.i_var = i_var;
+    else
+        out.plot_shock_decomposition_info.i_var = unique([i_var(:); out.plot_shock_decomposition_info.i_var(:)]);
+    end
+end
+
+type=options_.plot_shock_decomp.type;
+if isequal(type, 'aoa') && isfield(options_.plot_shock_decomp,'q2a') && isstruct(options_.plot_shock_decomp.q2a)
+    q2avec=options_.plot_shock_decomp.q2a;
+    if nvar>1
+        for jv = 1:nvar
+            my_varlist = varlist(jv);
+            indv = strcmp(my_varlist,{q2avec.qname});
+            options_.plot_shock_decomp.q2a =  q2avec(indv);
+            plot_shock_decomposition(M_,oo_,options_,my_varlist);
+        end
+        return
+    else
+        indv = strcmp(varlist,{q2avec.qname});
+        options_.plot_shock_decomp.q2a =  q2avec(indv);
+    end
+end
+
 % number of variables
 endo_nbr = M_.endo_nbr;
 
@@ -76,10 +146,14 @@ if ~isfield(options_.plot_shock_decomp,'init_cond_decomp')
     options_.plot_shock_decomp.init_cond_decomp=0;
 end
 
+options_.plot_shock_decomp.initval=0;
 if ~isempty(options_.plot_shock_decomp.fig_name)
     fig_name=[' ' options_.plot_shock_decomp.fig_name];
+    if length(fig_name)>=8 && strcmp(fig_name(end-6:end),'initval')
+        options_.plot_shock_decomp.initval=1;
+    end
 end
-type=options_.plot_shock_decomp.type;
+
 detail_plot=options_.plot_shock_decomp.detail_plot;
 realtime_= options_.plot_shock_decomp.realtime;
 vintage_ = options_.plot_shock_decomp.vintage;
@@ -200,8 +274,8 @@ if isequal(type,'aoa') && isstruct(q2a)
         if isempty(t0)
             error('the realtime decompositions are not stored in Q4! Please check your dates and settings.')
         end
-        if ~isfield(q2a,'var_type') % private trap for aoa calls
-            q2a.var_type=1;
+        if ~isfield(q2a,'type') % private trap for aoa calls
+            q2a.type=1;
         end
         if ~isfield(q2a,'islog') % private trap for aoa calls
             q2a.islog=0;
@@ -251,6 +325,20 @@ if options_.plot_shock_decomp.use_shock_groups
             [z, shock_names, M_] = make_the_groups(z,gend,endo_nbr,nshocks,M_,options_);
             M_.endo_names = endo_names;
             M_.endo_names_tex = endo_names_tex;
+        else
+            % here we know we only have one variable to handle
+            if isstruct(q2a.aux) && ischar(q2a.aux.y)
+                steady_state_aux  = get_mean(q2a.aux.y);
+                q2a.aux.y=repmat(steady_state_aux,16,1);
+                q2a.aux.yss=steady_state_aux;
+            end
+            [~, yssa, ~, gyssa] = ...
+                quarterly2annual(repmat(steady_state,16,1),steady_state,q2a.GYTREND0,q2a.type,q2a.islog,q2a.aux);
+            if q2a.plot==1
+                steady_state = gyssa;
+            else
+                steady_state = yssa;
+            end
         end
     else
         gend = size(z,3);
@@ -333,8 +421,8 @@ switch type
     end
     if isstruct(q2a)
         if realtime_ == 0
-            if ~isfield(q2a,'var_type') % private trap for aoa calls
-                q2a.var_type=1;
+            if ~isfield(q2a,'type') % private trap for aoa calls
+                q2a.type=1;
             end
             if ~isfield(q2a,'islog') % private trap for aoa calls
                 q2a.islog=0;
@@ -356,17 +444,19 @@ switch type
             if isstruct(q2a.aux) && ischar(q2a.aux.y)
                 opts=options_;
                 opts.plot_shock_decomp.type='qoq';
+                opts.plot_shock_decomp.use_shock_groups=[];
                 [y_aux, steady_state_aux] = plot_shock_decomposition(M_,oo_,opts,q2a.aux.y);
                 q2a.aux.y=y_aux;
                 q2a.aux.yss=steady_state_aux;
             end
+            i_var0 = i_var;
             [za, endo_names, endo_names_tex, steady_state, i_var, oo_] = ...
                 annualized_shock_decomposition(z,M_, options_, i_var, t0, options_.nobs, realtime_, vintage_, steady_state,q2a);
                 if options_.plot_shock_decomp.interactive && ~isempty(options_.plot_shock_decomp.use_shock_groups)
                     mygroup = options_.plot_shock_decomp.use_shock_groups;
                     options_.plot_shock_decomp.use_shock_groups='';
                     zafull = ...
-                        annualized_shock_decomposition(z,M_, options_, i_var, t0, options_.nobs, realtime_, vintage_, steady_state,q2a);
+                        annualized_shock_decomposition(zfull(i_var0,:,:),M_, options_, i_var, t0, options_.nobs, realtime_, vintage_, steady_state,q2a);
                     options_.plot_shock_decomp.use_shock_groups = mygroup;
                 end
             end
@@ -378,7 +468,7 @@ switch type
         M_.endo_names = endo_names;
         M_.endo_names_tex = endo_names_tex;
         %     endo_nbr = size(z,1);
-        if realtime_<2
+        if realtime_<2 || vintage_ == 0
             initial_date = initial_date1;
         else
             initial_date = initial_date0;
@@ -412,8 +502,8 @@ if steadystate
     options_.plot_shock_decomp.steady_state=steady_state;
 end
 
-if nargout
-    z=z(i_var,:,:);
+if nargout == 2
+    out=z(i_var,:,:);
     steady_state = steady_state(i_var);
     return
 end
@@ -427,7 +517,11 @@ if ~isempty(options_.plot_shock_decomp.plot_init_date)
     a = find((initial_date:initial_date+b-1)==options_.plot_shock_decomp.plot_init_date);
 end
 if ~isempty(options_.plot_shock_decomp.plot_end_date)
-    b = find((initial_date:initial_date+b-1)==options_.plot_shock_decomp.plot_end_date);
+    if options_.plot_shock_decomp.plot_end_date<=(max(initial_date:initial_date+b-1))
+        b = find((initial_date:initial_date+b-1)==options_.plot_shock_decomp.plot_end_date);
+    else
+        warning('You set plot_end_date larger than smoother size!!');
+    end
 end
 z = z(:,:,a:b);
 % end crop data
@@ -438,10 +532,12 @@ if options_.plot_shock_decomp.interactive && ~isempty(options_.plot_shock_decomp
     options_.plot_shock_decomp.zfull = zfull;
 end
 
-if detail_plot
-    graph_decomp_detail(z, shock_names, M_.endo_names, i_var, my_initial_date, M_, options_)
-else
-    graph_decomp(z, shock_names, M_.endo_names, i_var, my_initial_date, M_, options_);
+if ~options_.no_graph.plot_shock_decomposition
+    if detail_plot
+        graph_decomp_detail(z, shock_names, M_.endo_names, i_var, my_initial_date, M_, options_);
+    else
+        graph_decomp(z, shock_names, M_.endo_names, i_var, my_initial_date, M_, options_);
+    end
 end
 
 if write_xls
diff --git a/matlab/realtime_shock_decomposition.m b/matlab/realtime_shock_decomposition.m
index 4160c10fc..f7574ab87 100644
--- a/matlab/realtime_shock_decomposition.m
+++ b/matlab/realtime_shock_decomposition.m
@@ -39,6 +39,18 @@ function oo_ = realtime_shock_decomposition(M_,oo_,options_,varlist,bayestopt_,e
 % You should have received a copy of the GNU General Public License
 % along with Dynare.  If not, see <http://www.gnu.org/licenses/>.
 
+
+if isfield(oo_,'shock_decomposition_info') && isfield(oo_.shock_decomposition_info,'i_var')
+    if isfield (oo_,'realtime_conditional_shock_decomposition') ...
+            || isfield (oo_,'realtime_forecast_shock_decomposition') ...
+            || isfield (oo_,'realtime_shock_decomposition') ...
+            || isfield (oo_,'shock_decomposition') ...
+            || isfield (oo_,'conditional_shock_decomposition') ...
+            || isfield (oo_,'initval_decomposition')
+        error('realtime_shock_decomposition::squeezed shock decompositions are already stored in oo_')
+    end
+end
+
 % indices of endogenous variables
 if isempty(varlist)
     varlist = M_.endo_names(1:M_.orig_endo_nbr);
@@ -68,9 +80,10 @@ if isempty(parameter_set)
     end
 end
 
-presample = max(1,options_.presample);
+presample = max(1,options_.presample-1);
 if isfield(options_.shock_decomp,'presample')
-    presample = max(presample,options_.shock_decomp.presample);
+    my_presample = max(1,options_.shock_decomp.presample);
+    presample = min(presample,my_presample);
 end
 % forecast_=0;
 forecast_ = options_.shock_decomp.forecast;
@@ -244,7 +257,9 @@ for j=presample+1:nobs
             %             zn(:,1:nshocks,i) = zn(:,1:nshocks,i) + B(inv_order_var,:).*repmat(epsilon(:,i+gend-forecast_-1)',endo_nbr,1);
             zn(:,nshocks+1,i) = zn(:,nshocks+2,i) - sum(zn(:,1:nshocks,i),2);
         end
-        oo_.conditional_shock_decomposition.(['time_' int2str(j-forecast_)])=zn;
+        if ismember(j-forecast_,save_realtime)
+            oo_.conditional_shock_decomposition.(['time_' int2str(j-forecast_)])=zn;
+        end
     end
     %%
 
@@ -260,28 +275,39 @@ for j=presample+1:nobs
 
     if forecast_
         zfrcst(:,:,j+1) = z(:,:,gend+1);
-        oo_.realtime_forecast_shock_decomposition.(['time_' int2str(j)])=z(:,:,gend:end);
+        ootmp.realtime_forecast_shock_decomposition.(['time_' int2str(j)])=z(:,:,gend:end);
+        if ismember(j,save_realtime)
+            oo_.realtime_forecast_shock_decomposition.(['time_' int2str(j)]) = ...
+                ootmp.realtime_forecast_shock_decomposition.(['time_' int2str(j)]);
+        end
         if j>forecast_+presample
             %% realtime conditional shock decomp k step ahead
-            oo_.realtime_conditional_shock_decomposition.(['time_' int2str(j-forecast_)]) = ...
+            ootmp.realtime_conditional_shock_decomposition.(['time_' int2str(j-forecast_)]) = ...
                 zreal(:,:,j-forecast_:j) - ...
-                oo_.realtime_forecast_shock_decomposition.(['time_' int2str(j-forecast_)]);
-            oo_.realtime_conditional_shock_decomposition.(['time_' int2str(j-forecast_)])(:,end-1,:) = ...
-                oo_.realtime_conditional_shock_decomposition.(['time_' int2str(j-forecast_)])(:,end-1,:) + ...
-                oo_.realtime_forecast_shock_decomposition.(['time_' int2str(j-forecast_)])(:,end,:);
-            oo_.realtime_conditional_shock_decomposition.(['time_' int2str(j-forecast_)])(:,end,:) = ...
+                ootmp.realtime_forecast_shock_decomposition.(['time_' int2str(j-forecast_)]);
+            ootmp.realtime_conditional_shock_decomposition.(['time_' int2str(j-forecast_)])(:,end-1,:) = ...
+                ootmp.realtime_conditional_shock_decomposition.(['time_' int2str(j-forecast_)])(:,end-1,:) + ...
+                ootmp.realtime_forecast_shock_decomposition.(['time_' int2str(j-forecast_)])(:,end,:);
+            ootmp.realtime_conditional_shock_decomposition.(['time_' int2str(j-forecast_)])(:,end,:) = ...
                 zreal(:,end,j-forecast_:j);
-
+            if ismember(j-forecast_,save_realtime)
+                oo_.realtime_conditional_shock_decomposition.(['time_' int2str(j-forecast_)]) = ...
+                    ootmp.realtime_conditional_shock_decomposition.(['time_' int2str(j-forecast_)]);
+            end
             if j==nobs
                 for my_forecast_=(forecast_-1):-1:1
-                    oo_.realtime_conditional_shock_decomposition.(['time_' int2str(j-my_forecast_)]) = ...
+                    ootmp.realtime_conditional_shock_decomposition.(['time_' int2str(j-my_forecast_)]) = ...
                         zreal(:,:,j-my_forecast_:j) - ...
-                        oo_.realtime_forecast_shock_decomposition.(['time_' int2str(j-my_forecast_)])(:,:,1:my_forecast_+1);
-                    oo_.realtime_conditional_shock_decomposition.(['time_' int2str(j-my_forecast_)])(:,end-1,:) = ...
-                        oo_.realtime_conditional_shock_decomposition.(['time_' int2str(j-my_forecast_)])(:,end-1,:) + ...
-                        oo_.realtime_forecast_shock_decomposition.(['time_' int2str(j-my_forecast_)])(:,end,1:my_forecast_+1);
-                    oo_.realtime_conditional_shock_decomposition.(['time_' int2str(j-my_forecast_)])(:,end,:) = ...
+                        ootmp.realtime_forecast_shock_decomposition.(['time_' int2str(j-my_forecast_)])(:,:,1:my_forecast_+1);
+                    ootmp.realtime_conditional_shock_decomposition.(['time_' int2str(j-my_forecast_)])(:,end-1,:) = ...
+                        ootmp.realtime_conditional_shock_decomposition.(['time_' int2str(j-my_forecast_)])(:,end-1,:) + ...
+                        ootmp.realtime_forecast_shock_decomposition.(['time_' int2str(j-my_forecast_)])(:,end,1:my_forecast_+1);
+                    ootmp.realtime_conditional_shock_decomposition.(['time_' int2str(j-my_forecast_)])(:,end,:) = ...
                         zreal(:,end,j-my_forecast_:j);
+                    if ismember(j-my_forecast_,save_realtime)
+                        oo_.realtime_conditional_shock_decomposition.(['time_' int2str(j-my_forecast_)]) = ...
+                            ootmp.realtime_conditional_shock_decomposition.(['time_' int2str(j-my_forecast_)]);
+                    end
                 end
             end
 
diff --git a/matlab/set_default_option.m b/matlab/set_default_option.m
index 43910ba29..a388ca94e 100644
--- a/matlab/set_default_option.m
+++ b/matlab/set_default_option.m
@@ -41,7 +41,7 @@ if isempty(options.(field))
     return
 end
 
-if ~iscell(options.(field)) && ~isdates(options.(field))
+if ~iscell(options.(field)) && ~isdates(options.(field)) && ~isstruct(options.(field))
     if isnan(options.(field))
         options.(field) = default;
         return
diff --git a/matlab/shock_decomposition.m b/matlab/shock_decomposition.m
index 2793b1a9e..20e6adfe8 100644
--- a/matlab/shock_decomposition.m
+++ b/matlab/shock_decomposition.m
@@ -41,6 +41,16 @@ function [oo_,M_] = shock_decomposition(M_,oo_,options_,varlist,bayestopt_,estim
 % along with Dynare.  If not, see <http://www.gnu.org/licenses/>.
 
 % indices of endogenous variables
+
+if isfield(oo_,'shock_decomposition_info') && isfield(oo_.shock_decomposition_info,'i_var')
+    if isfield (oo_,'realtime_conditional_shock_decomposition') ...
+            || isfield (oo_,'realtime_forecast_shock_decomposition') ...
+            || isfield (oo_,'realtime_shock_decomposition') ...
+            || isfield (oo_,'conditional_shock_decomposition') ...
+            || isfield (oo_,'initval_decomposition')
+        error('shock_decomposition::squeezed shock decompositions are already stored in oo_')
+    end
+end
 if isempty(varlist)
     varlist = M_.endo_names(1:M_.orig_endo_nbr);
 end
@@ -123,5 +133,5 @@ end
 oo_.shock_decomposition = z;
 
 if ~options_.no_graph.shock_decomposition
-    plot_shock_decomposition(M_,oo_,options_,varlist);
+    oo_ = plot_shock_decomposition(M_,oo_,options_,varlist);
 end
\ No newline at end of file
diff --git a/matlab/squeeze_shock_decomp.m b/matlab/squeeze_shock_decomp.m
new file mode 100644
index 000000000..49073ee08
--- /dev/null
+++ b/matlab/squeeze_shock_decomp.m
@@ -0,0 +1,63 @@
+function oo_ = squeeze_shock_decomp(M_,oo_,options_,var_list_)
+
+if isfield(oo_,'plot_shock_decomposition_info') && isfield(oo_.plot_shock_decomposition_info','i_var')
+    my_vars = oo_.plot_shock_decomposition_info.i_var;
+else
+    my_vars=[];
+end
+if nargin>3
+    my_vars = [varlist_indices(var_list_,M_.endo_names); my_vars];
+end
+sd_vlist = M_.endo_names(my_vars,:);
+
+if isfield(options_.plot_shock_decomp,'q2a') && isstruct(options_.plot_shock_decomp.q2a)
+    
+    avname={options_.plot_shock_decomp.q2a.qname};
+    sda = options_.plot_shock_decomp.q2a(ismember(avname,sd_vlist));
+    for k=1:length(sda)
+        if isstruct(sda(k).aux)
+            sd_vlist = [sd_vlist; {sda(k).aux.y}];
+        end
+    end
+end
+
+if isempty(sd_vlist)
+    disp('Nothing has been squeezed: there is no list of variables for it!')
+    return
+end
+i_var = varlist_indices(sd_vlist,M_.endo_names);
+
+oo_.plot_shock_decomposition_info.i_var = i_var;
+oo_.shock_decomposition_info.i_var = i_var;
+if isfield (oo_,'shock_decomposition')
+    oo_.shock_decomposition = oo_.shock_decomposition(i_var,:,:);
+end
+if isfield (oo_,'realtime_conditional_shock_decomposition')
+    oo_.realtime_conditional_shock_decomposition = ...
+        my_squeeze(oo_.realtime_conditional_shock_decomposition, i_var);
+end
+if isfield (oo_,'realtime_forecast_shock_decomposition')
+    oo_.realtime_forecast_shock_decomposition = ...
+        my_squeeze(oo_.realtime_forecast_shock_decomposition, i_var);
+end
+if isfield (oo_,'realtime_shock_decomposition')
+    oo_.realtime_shock_decomposition = ...
+        my_squeeze(oo_.realtime_shock_decomposition, i_var);
+end
+if isfield (oo_,'conditional_shock_decomposition')
+    oo_.conditional_shock_decomposition = ...
+        my_squeeze(oo_.conditional_shock_decomposition, i_var);
+end
+if isfield (oo_,'initval_decomposition')
+    oo_.initval_decomposition = oo_.initval_decomposition(i_var,:,:);
+end
+
+end
+
+function shock_decomposition = my_squeeze(shock_decomposition, i_var)
+fnam = fieldnames(shock_decomposition);
+for k=1:length(fnam)
+    shock_decomposition.(fnam{k}) =  shock_decomposition.(fnam{k})(i_var,:,:);
+end
+
+end
diff --git a/tests/shock_decomposition/ls2003_plot.mod b/tests/shock_decomposition/ls2003_plot.mod
index e07876b20..3ef63d277 100644
--- a/tests/shock_decomposition/ls2003_plot.mod
+++ b/tests/shock_decomposition/ls2003_plot.mod
@@ -83,10 +83,14 @@ A e_A;
 end;
 
 options_.initial_date=dates('1989Q4'); % date arbitrarily set for testing purposes
-shock_decomposition(use_shock_groups=trade) y_obs R_obs pie_obs dq de;
+shock_decomposition(nograph);
+// test for nothing to squeeze
+oo_ = squeeze_shock_decomp(M_,oo_,options_);
 
 // standard plot
 plot_shock_decomposition y_obs R_obs pie_obs dq de;
+// grouped shocks
+plot_shock_decomposition(use_shock_groups=trade) y_obs R_obs pie_obs dq de;
 
 // test datailed, custom name and yoy plots
 plot_shock_decomposition(detail_plot, fig_name = MR, type = yoy) y_obs R_obs pie_obs dq de;
@@ -107,6 +111,9 @@ close all,
 // first compute realtime decompositions [pre-processor not yet available]
 realtime_shock_decomposition(forecast=8, save_realtime=[5 9 13 17 21 25 29 33 37 41 45 49 53 57 61 65 69 73 77]) y_obs R_obs pie_obs dq de;
 
+// test squeeze
+oo_ = squeeze_shock_decomp(M_,oo_,options_);
+
 //realtime pooled
 plot_shock_decomposition(realtime = 1) y_obs R_obs pie_obs dq de;
 
@@ -129,6 +136,7 @@ close all,
 // now I test annualized variables
 // options_.plot_shock_decomp.q2a=1;
 // options_.plot_shock_decomp.islog=1;
+// this also triggers re-computing of decompositions since y was not present in squeeze set
 plot_shock_decomposition(detail_plot, type = aoa) y;
 
 plot_shock_decomposition(realtime = 1) y;
@@ -148,6 +156,9 @@ close all,
 // testing realtime decomposition with fast_realtime option
 realtime_shock_decomposition(fast_realtime=75) y_obs R_obs pie_obs dq de;
 
+// re-test squeeze
+oo_ = squeeze_shock_decomp(M_,oo_,options_);
+
 collect_latex_files;
 if system(['pdflatex -halt-on-error -interaction=batchmode ' M_.fname '_TeX_binder.tex'])
     error('TeX-File did not compile.')