1function [zout, zss] = epilogue_shock_decomposition(zout ,M_, oo_)
2
3ivar = varlist_indices(M_.epilogue_var_list_,M_.endo_names);
4y = nan(length(M_.epilogue_names),size(zout,2),size(zout,3));
5z=zout(ivar,:,:);
6if isfield(oo_.dr,'ys')
7    zss = oo_.dr.ys;
8else
9    zss = oo_.steady_state;
10end
11
12ztmp = dseries(zss(ivar),'',M_.epilogue_var_list_);
13fname = M_.fname;
14h_dynamic = str2func([fname '.epilogue_dynamic']);
15h_static = str2func([fname '.epilogue_static']);
16yss = extract(h_static(M_.params, ztmp),M_.epilogue_names{:});
17yss = yss.data;
18yss1 = repmat(yss(:),[1,1,size(y,3)]);
19for k=1:size(z,2)
20    ztmp = dseries(squeeze(z(:,k,:))+zss(ivar),'',M_.epilogue_var_list_);
21    tmp = extract(h_dynamic(M_.params, ztmp),M_.epilogue_names{:});
22    y(:,k,:) = tmp.data';
23    y(:,k,:) = y(:,k,:)-yss1;
24end
25
26nterms = size(z,2);
27for k=1:size(y,1)
28    ytmp  = squeeze(y(k,:,:));
29    yres = ytmp(end,:) - sum(ytmp(1:end-1,:));
30    if ~isoctave && matlab_ver_less_than('9.1') % Automatic broadcasting was introduced in MATLAB R2016b
31        w = bsxfun(@rdivide, abs(ytmp(1:end-1,:)), sum(abs(ytmp(1:end-1,:))))
32    else
33        w = abs(ytmp(1:end-1,:))./sum(abs(ytmp(1:end-1,:)));
34    end
35    %     ytmp(1:end-1,:) = ytmp(1:end-1,:) + repmat(yres,[nterms-1 1])/(nterms-1);
36    ytmp(1:end-1,:) = ytmp(1:end-1,:) + repmat(yres,[nterms-1 1]).*w;
37    y(k,:,:) = ytmp;
38end
39
40zout = cat(1,zout,y);
41zss = [zss; yss(:)];
42