diff --git a/matlab/annualized_shock_decomposition.m b/matlab/annualized_shock_decomposition.m
index 03b5feb443b9ce6d083c59c18fb4e30d5be10735..038f15a684f7a9401000a40b65bc7876605bcfbd 100644
--- a/matlab/annualized_shock_decomposition.m
+++ b/matlab/annualized_shock_decomposition.m
@@ -305,16 +305,22 @@ for j=1:nvar
     end
     ztmp=squeeze(za(j,:,:));
     if cumfix==0
-        zscale = sum(ztmp(1:end-1,:))./ztmp(end,:);
-        ztmp(1:end-1,:) = ztmp(1:end-1,:)./repmat(zscale,[nterms-1,1]);
+        zres = ztmp(end,:) - sum(ztmp(1:end-1,:));
+        w = abs(ztmp(1:end-1,:))./sum(abs(ztmp(1:end-1,:)));
+        ztmp(1:end-1,:) = ztmp(1:end-1,:) + repmat(zres,[nterms-1 1]).*w;
+%         zscale = sum(ztmp(1:end-1,:))./ztmp(end,:);
+%         ztmp(1:end-1,:) = ztmp(1:end-1,:)./repmat(zscale,[nterms-1,1]);
     else
         zres = ztmp(end,:)-sum(ztmp(1:end-1,:));
         ztmp(1:end-1,:) = ztmp(1:end-1,:) + repmat(zres,[nterms-1 1])/(nterms-1);
     end
     gztmp=squeeze(gza(j,:,:));
     if cumfix==0
-        gscale = sum(gztmp(1:end-1,:))./ gztmp(end,:);
-        gztmp(1:end-1,:) = gztmp(1:end-1,:)./repmat(gscale,[nterms-1,1]);
+        gres = gztmp(end,:) - sum(gztmp(1:end-1,:));
+        w = abs(gztmp(1:end-1,:))./sum(abs(gztmp(1:end-1,:)));
+        gztmp(1:end-1,:) = gztmp(1:end-1,:) + repmat(gres,[nterms-1 1]).*w;
+%         gscale = sum(gztmp(1:end-1,:))./ gztmp(end,:);
+%         gztmp(1:end-1,:) = gztmp(1:end-1,:)./repmat(gscale,[nterms-1,1]);
     else
         gres = gztmp(end,:) - sum(gztmp(1:end-1,:));
         gztmp(1:end-1,:) = gztmp(1:end-1,:) + repmat(gres,[nterms-1 1])/(nterms-1);
diff --git a/matlab/default_option_values.m b/matlab/default_option_values.m
index d85cee7e6a003c50b199fb135c1df9546c00eedc..c64d246e851bfd6c2a4d177c9348af081b744d98 100644
--- a/matlab/default_option_values.m
+++ b/matlab/default_option_values.m
@@ -652,6 +652,7 @@ options_.parameter_set = [];
 options_.use_shock_groups = '';
 options_.shock_decomp.colormap = '';
 options_.shock_decomp.init_state = 0;
+options_.shock_decomp.with_epilogue = false;
 
 % Shock decomposition realtime
 options_.shock_decomp.forecast = 0;
diff --git a/matlab/epilogue_shock_decomposition.m b/matlab/epilogue_shock_decomposition.m
new file mode 100644
index 0000000000000000000000000000000000000000..51f5cbc406ea8079dbd3ec2f2da674172160f0df
--- /dev/null
+++ b/matlab/epilogue_shock_decomposition.m
@@ -0,0 +1,39 @@
+function [zout, zss] = epilogue_shock_decomposition(zout ,M_, oo_)
+
+ivar = varlist_indices(M_.epilogue_var_list_,M_.endo_names);
+y = nan(length(M_.epilogue_names),size(zout,2),size(zout,3));
+z=zout(ivar,:,:);
+if isfield(oo_.dr,'ys')
+    zss = oo_.dr.ys;
+else
+    zss = oo_.steady_state;
+end
+
+ztmp = dseries(zss(ivar),'',M_.epilogue_var_list_);
+fname = M_.fname;
+h_dynamic = str2func([fname '.epilogue_dynamic']);
+h_static = str2func([fname '.epilogue_static']);
+yss = extract(h_static(M_.params, ztmp),M_.epilogue_names{:});
+yss = yss.data;
+yss1 = repmat(yss(:),[1,1,size(y,3)]);
+for k=1:size(z,2)
+    ztmp = dseries(squeeze(z(:,k,:))+zss(ivar),'',M_.epilogue_var_list_);
+    tmp = extract(h_dynamic(M_.params, ztmp),M_.epilogue_names{:});
+    y(:,k,:) = tmp.data';
+    y(:,k,:) = y(:,k,:)-yss1;
+end
+
+nterms = size(z,2);
+for k=1:size(y,1)
+    ytmp  = squeeze(y(k,:,:));
+    yres = ytmp(end,:) - sum(ytmp(1:end-1,:));
+    w = abs(ytmp(1:end-1,:))./sum(abs(ytmp(1:end-1,:)));
+%     ytmp(1:end-1,:) = ytmp(1:end-1,:) + repmat(yres,[nterms-1 1])/(nterms-1);
+    ytmp(1:end-1,:) = ytmp(1:end-1,:) + repmat(yres,[nterms-1 1]).*w;
+    y(k,:,:) = ytmp;
+end
+
+zout = cat(1,zout,y);
+zss = [zss; yss(:)];
+
+
diff --git a/matlab/initial_condition_decomposition.m b/matlab/initial_condition_decomposition.m
index db35777ae18d7460fe5171dcfdf12f82a5eff3ac..9ec20eda6c478b8005fc4df8a28f3cce36115406 100644
--- a/matlab/initial_condition_decomposition.m
+++ b/matlab/initial_condition_decomposition.m
@@ -96,6 +96,7 @@ if ~isfield(oo_,'initval_decomposition') || isequal(varlist,0)
             error('initval_decomposition::squeezed shock decompositions are already stored in oo_')
         end
     end
+    with_epilogue = options_.initial_condition_decomp.with_epilogue;
     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_);
@@ -141,6 +142,12 @@ if ~isfield(oo_,'initval_decomposition') || isequal(varlist,0)
 
     end
 
+    if with_epilogue
+        [z, epilogue_steady_state] = epilogue_shock_decomposition(z, M_, oo_);
+        if ~isfield(oo_,'shock_decomposition_info') || ~isfield(oo_.shock_decomposition_info,'epilogue_steady_state')
+            oo_.shock_decomposition_info.epilogue_steady_state = epilogue_steady_state;
+        end
+    end
     oo_.initval_decomposition = z;
 end
 
diff --git a/matlab/plot_shock_decomposition.m b/matlab/plot_shock_decomposition.m
index 50b978cbcd9bd9023213af4542761b73d1c235d4..5c35fa2246f21021cac8ebf43715254d0bba19f7 100644
--- a/matlab/plot_shock_decomposition.m
+++ b/matlab/plot_shock_decomposition.m
@@ -59,11 +59,26 @@ if ~isempty(init2shocks)
     init2shocks=M_.init2shocks.(init2shocks);
 end
 
+
+epilogue_decomp=false;
+if exist_varlist && any(ismember(varlist,M_.epilogue_names))
+    epilogue_decomp=true;
+    M_.endo_names = [M_.endo_names;M_.epilogue_names];
+    M_.endo_names_tex = [M_.endo_names_tex;M_.epilogue_names];
+    M_.endo_nbr = length( M_.endo_names );
+end
 if isfield(oo_.shock_decomposition_info,'i_var') && (M_.endo_nbr>=M_.orig_endo_nbr)
+    if max(oo_.shock_decomposition_info.i_var)>M_.orig_endo_nbr
+        epilogue_decomp=true;
+        M_.endo_names = [M_.endo_names;M_.epilogue_names];
+        M_.endo_names_tex = [M_.endo_names_tex;M_.epilogue_names];
+        M_.endo_nbr = length( M_.endo_names );
+    end
     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
@@ -119,9 +134,6 @@ if isequal(type, 'aoa') && isfield(options_.plot_shock_decomp,'q2a') && isstruct
     end
 end
 
-% number of variables
-endo_nbr = M_.endo_nbr;
-
 % number of shocks
 nshocks = M_.exo_nbr;
 fig_name='';
@@ -253,10 +265,14 @@ if ~isempty(init2shocks) && ~expand
     end   
 end    
 
-if isfield(oo_.dr,'ys')
-    steady_state = oo_.dr.ys;
+if ~epilogue_decomp
+    if isfield(oo_.dr,'ys')
+        steady_state = oo_.dr.ys;
+    else
+        steady_state = oo_.steady_state;
+    end
 else
-    steady_state = oo_.steady_state;
+    steady_state = oo_.shock_decomposition_info.epilogue_steady_state;
 end
 
 if isequal(type,'aoa') && isstruct(q2a)
@@ -343,6 +359,7 @@ if options_.plot_shock_decomp.use_shock_groups
     else
         gend = size(z,3);
         zfull = z;
+        endo_nbr = size(z,1);
         [z, shock_names, M_] = make_the_groups(z,gend,endo_nbr,nshocks,M_,options_);
     end
     if ~isempty(init2shocks) && ~expand
diff --git a/matlab/realtime_shock_decomposition.m b/matlab/realtime_shock_decomposition.m
index f7574ab8759f036b3182fbe2971f313935f43f64..e003710da0122bb8382ebfd9c812a9b002a1207b 100644
--- a/matlab/realtime_shock_decomposition.m
+++ b/matlab/realtime_shock_decomposition.m
@@ -51,6 +51,8 @@ if isfield(oo_,'shock_decomposition_info') && isfield(oo_.shock_decomposition_in
     end
 end
 
+with_epilogue = options_.shock_decomp.with_epilogue;
+
 % indices of endogenous variables
 if isempty(varlist)
     varlist = M_.endo_names(1:M_.orig_endo_nbr);
@@ -100,8 +102,8 @@ end
 save_realtime = options_.shock_decomp.save_realtime;
 % array of time points in the range options_.presample+1:options_.nobs
 
-zreal = zeros(endo_nbr,nshocks+2,options_.nobs+forecast_);
-zcond = zeros(endo_nbr,nshocks+2,options_.nobs);
+zreal = zeros(endo_nbr+length(M_.epilogue_names)*with_epilogue,nshocks+2,options_.nobs+forecast_);
+zcond = zeros(endo_nbr+length(M_.epilogue_names)*with_epilogue,nshocks+2,options_.nobs);
 
 options_.selected_variables_only = 0; %make sure all variables are stored
 options_.plot_priors=0;
@@ -226,6 +228,12 @@ for j=presample+1:nobs
         z(:,nshocks+1,i) = z(:,nshocks+2,i) - sum(z(:,1:nshocks,i),2);
     end
 
+    if with_epilogue
+        [z, epilogue_steady_state] = epilogue_shock_decomposition(z, M_, oo_);
+        if ~isfield(oo_,'shock_decomposition_info') || ~isfield(oo_.shock_decomposition_info,'epilogue_steady_state')
+            oo_.shock_decomposition_info.epilogue_steady_state = epilogue_steady_state;
+        end
+    end
     %% conditional shock decomp 1 step ahead
     z1 = zeros(endo_nbr,nshocks+2);
     z1(:,end) = Smoothed_Variables_deviation_from_mean(:,gend);
@@ -234,6 +242,15 @@ for j=presample+1:nobs
         z1(:,1:nshocks) = z1(:,1:nshocks) + B(inv_order_var,:).*repmat(epsilon(:,i)',endo_nbr,1);
         z1(:,nshocks+1) = z1(:,nshocks+2) - sum(z1(:,1:nshocks),2);
     end
+    if with_epilogue
+        clear ztmp0
+        ztmp0(:,1,:) = Smoothed_Variables_deviation_from_mean(:,1:gend-1);
+        ztmp0(:,2,:) = Smoothed_Variables_deviation_from_mean(:,1:gend-1);
+        ztmp = cat(3,cat(2,zeros(endo_nbr,nshocks,gend-1),ztmp0),z1);
+%         ztmp = cat(3,zeros(endo_nbr,nshocks+2,40),ztmp); % pad with zeros in presample
+        z1  = epilogue_shock_decomposition(ztmp, M_, oo_);
+        z1=squeeze(z1(:,:,end));
+    end
     %%
 
     %% conditional shock decomp k step ahead
@@ -257,6 +274,15 @@ 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
+        if with_epilogue
+            clear ztmp0
+            ztmp0(:,1,:) = Smoothed_Variables_deviation_from_mean(:,1:gend-forecast_-1);
+            ztmp0(:,2,:) = Smoothed_Variables_deviation_from_mean(:,1:gend-forecast_-1);
+            ztmp = cat(3,cat(2,zeros(endo_nbr,nshocks,gend-forecast_-1),ztmp0),zn);
+%             ztmp = cat(3,zeros(endo_nbr,nshocks+2,40),ztmp); % pad with zeros (st state) in presample
+            zn  = epilogue_shock_decomposition(ztmp, M_, oo_);
+            zn=squeeze(zn(:,:,end-forecast_:end));
+        end
         if ismember(j-forecast_,save_realtime)
             oo_.conditional_shock_decomposition.(['time_' int2str(j-forecast_)])=zn;
         end
@@ -265,8 +291,10 @@ for j=presample+1:nobs
 
     if init
         zreal(:,:,1:j) = z(:,:,1:j);
-    else
+    elseif j<nobs
         zreal(:,:,j) = z(:,:,gend);
+    else
+        zreal(:,:,j:end) = z(:,:,gend:end);
     end
     zcond(:,:,j) = z1;
     if ismember(j,save_realtime)
diff --git a/matlab/set_default_initial_condition_decomposition_options.m b/matlab/set_default_initial_condition_decomposition_options.m
index 167addad5a5c3bbad37484a6ba0e3f7ecdc9da2a..f354db07da953ab365c2ddaa85be9c3306fe773d 100644
--- a/matlab/set_default_initial_condition_decomposition_options.m
+++ b/matlab/set_default_initial_condition_decomposition_options.m
@@ -35,6 +35,7 @@ options.initial_condition_decomp.fig_name = '';
 options.initial_condition_decomp.detail_plot = false;
 options.initial_condition_decomp.init2shocks = [];
 options.initial_condition_decomp.steadystate = false;
+options.initial_condition_decomp.with_epilogue = options.shock_decomp.with_epilogue;
 options.initial_condition_decomp.write_xls = false;
 options.initial_condition_decomp.type = '';
 options.initial_condition_decomp.plot_init_date = [];
diff --git a/matlab/shock_decomposition.m b/matlab/shock_decomposition.m
index 20e6adfe8596901898812687a89eef0aa6ccf913..274e5933485357cf4993829b625f92a85e642838 100644
--- a/matlab/shock_decomposition.m
+++ b/matlab/shock_decomposition.m
@@ -51,6 +51,8 @@ if isfield(oo_,'shock_decomposition_info') && isfield(oo_.shock_decomposition_in
         error('shock_decomposition::squeezed shock decompositions are already stored in oo_')
     end
 end
+with_epilogue = options_.shock_decomp.with_epilogue;
+
 if isempty(varlist)
     varlist = M_.endo_names(1:M_.orig_endo_nbr);
 end
@@ -130,8 +132,12 @@ for i=1:gend
     z(:,nshocks+1,i) = z(:,nshocks+2,i) - sum(z(:,1:nshocks,i),2);
 end
 
+if with_epilogue
+    [z, oo_.shock_decomposition_info.epilogue_steady_state] = epilogue_shock_decomposition(z, M_, oo_);
+end
+
 oo_.shock_decomposition = z;
 
 if ~options_.no_graph.shock_decomposition
     oo_ = plot_shock_decomposition(M_,oo_,options_,varlist);
-end
\ No newline at end of file
+end
diff --git a/matlab/squeeze_shock_decomposition.m b/matlab/squeeze_shock_decomposition.m
index ff01875a2d89a79c26c365e179d54abf7056fc50..56513c7174edfbf0761b42edc4a2aeba663c8399 100644
--- a/matlab/squeeze_shock_decomposition.m
+++ b/matlab/squeeze_shock_decomposition.m
@@ -29,15 +29,21 @@ function oo_ = squeeze_shock_decomposition(M_,oo_,options_,var_list_)
 % You should have received a copy of the GNU General Public License
 % along with Dynare.  If not, see <http://www.gnu.org/licenses/>.
 
+if ~options_.shock_decomp.with_epilogue
+    endo_names = M_.endo_names;
+else
+    endo_names = [M_.endo_names; M_.epilogue_names];
+    
+end
 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];
+    my_vars = [varlist_indices(var_list_,endo_names); my_vars];
 end
-sd_vlist = M_.endo_names(my_vars,:);
+sd_vlist = endo_names(my_vars,:);
 
 if isfield(options_.plot_shock_decomp,'q2a') && isstruct(options_.plot_shock_decomp.q2a)
     
@@ -54,7 +60,7 @@ 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);
+i_var = varlist_indices(sd_vlist,endo_names);
 
 oo_.plot_shock_decomposition_info.i_var = i_var;
 oo_.shock_decomposition_info.i_var = i_var;
diff --git a/tests/shock_decomposition/ls2003_plot.mod b/tests/shock_decomposition/ls2003_plot.mod
index a3134d2faad37e8e9ac4e2d3ceac2aaa1cd121a0..af0a05ff80671003f702a33b774ba748bad109bb 100644
--- a/tests/shock_decomposition/ls2003_plot.mod
+++ b/tests/shock_decomposition/ls2003_plot.mod
@@ -31,6 +31,13 @@ pie_obs = 4*pie;
 R_obs = 4*R;
 end;
 
+epilogue;
+// annualized level of y
+ya = exp(y)+exp(y(-1))+exp(y(-2))+exp(y(-3));
+// annualized growth rate of y
+gya = ya/ya(-4)-1;
+end;
+
 shocks;
 var e_R = 1.25^2;
 var e_q = 2.5^2;
@@ -83,9 +90,11 @@ A e_A;
 end;
 
 options_.initial_date=dates('1989Q4'); % date arbitrarily set for testing purposes
+options_.shock_decomp.with_epilogue=true;
 shock_decomposition(nograph);
 // test for nothing to squeeze
 squeeze_shock_decomposition;
+// oo_ = squeeze_shock_decomp(M_,oo_,options_);
 
 // standard plot
 plot_shock_decomposition y_obs R_obs pie_obs dq de;
@@ -113,6 +122,7 @@ realtime_shock_decomposition(forecast=8, save_realtime=[5 9 13 17 21 25 29 33 37
 
 // test squeeze
 squeeze_shock_decomposition;
+// oo_ = squeeze_shock_decomp(M_,oo_,options_);
 
 //realtime pooled
 plot_shock_decomposition(realtime = 1) y_obs R_obs pie_obs dq de;
@@ -134,10 +144,9 @@ plot_shock_decomposition(detail_plot, realtime = 3, vintage = 29) y_obs R_obs pi
 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
+    // 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(detail_plot, type = aoa) ya;
 
 plot_shock_decomposition(realtime = 1) y;
 plot_shock_decomposition(realtime = 1, vintage = 29) y;
@@ -158,6 +167,7 @@ realtime_shock_decomposition(fast_realtime=75) y_obs R_obs pie_obs dq de;
 
 // re-test squeeze
 squeeze_shock_decomposition;
+// oo_ = squeeze_shock_decomp(M_,oo_,options_);
 
 collect_latex_files;
 if system(['pdflatex -halt-on-error -interaction=batchmode ' M_.fname '_TeX_binder.tex'])