% Outputs steady state or HP-filtered simulated variables
function [SS IRF SIM] = dyn_extract(M_,oo_,oo_sim)

    % Extract steady state as structure
    SS = cell2struct(...
      num2cell(oo_.steady_state),...
      cellstr(M_.endo_names)...
      );

    % Extract impulse responses as structure
    IRF = oo_.irfs;
    
    if nargout>2

      % Remove first obeservations from simulation
      [xxx, V N] = size(oo_sim); %#ok<ASGLU>
      remove = 0; 
      oo_sim = reshape(oo_sim(:,:,remove+1:N),V,N-remove);

      % Convert to structure
      SIM = cell2struct(num2cell(oo_sim,2),cellstr(M_.endo_names));

      % Real variables;
      SIM.rexH  = SIM.xF*SS.PH;
      SIM.rexF  = SIM.xH*SS.PF;
      SIM.rgdpH = SIM.cH + SIM.iH + SIM.rexH - SIM.rexF;
      SIM.rgdpF = SIM.cF + SIM.iF + SIM.rexF - SIM.rexH;
      SIM.rnxH  = (SIM.rexH - SIM.rexF)./SIM.rgdpH;
      SIM.rnxF  = (SIM.rexF - SIM.rexH)./SIM.rgdpF;

      % Apply HP filter
      SIM = structfun(@(x) myhpfilter(x),SIM,'UniformOutput',false);
    
    end
    
end

% Takes log if all variables are nonnegative,
% then extracts cyclical component of HP filter
function w = myhpfilter(x)
  if all(x>=0), x = log(x); end
  diag = eye(length(x));
  d    = diff(diag,2); 
  w    = x' - (diag + 1600*(d'*d))\x'; 
%  w = x' - hpfilter(x,1600);
end