diff --git a/matlab/default_option_values.m b/matlab/default_option_values.m
index c4512438316d5e3ba34b387db5d4d179ab7cddd8..cede9d669827782ffb32b22cb51d7e7aa8e90536 100644
--- a/matlab/default_option_values.m
+++ b/matlab/default_option_values.m
@@ -651,6 +651,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 3fbfa13617c6b4fadd55bb776bed339472ba367d..63cdcbef2aaaf0ba4fa9fa44b5a1d62b7c0e1e8c 100644
--- a/matlab/squeeze_shock_decomposition.m
+++ b/matlab/squeeze_shock_decomposition.m
@@ -1,14 +1,20 @@
 function oo_ = squeeze_shock_decomposition(M_,oo_,options_,var_list_)
 
+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)
     
@@ -25,7 +31,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;