v4.1: added functions for shock_decomposition consistent with smoother

Dynare interface is missing
      need to add mechanism for aggregating shocks


git-svn-id: https://www.dynare.org/svn/dynare/trunk@2525 ac1d8469-bf42-47a9-8791-bf33cf982152

matlab/graph_decomp.m

@@ -0,0 +1,60 @@
+function []=graph_decomp(z,varlist,initial_period,freq)
+    global M_ 
+    exo_nbr = M_.exo_nbr;
+
+    gend = size(z,3);
+    x = initial_period-1/freq:(1/freq):initial_period+(gend-1)/freq;
+    
+    [i_var,nvar] = varlist_indices(varlist);
+    for j=1:nvar;
+        z1 = squeeze(z(i_var(j),:,:));
+        xmin = x(1);
+        xmax = x(end);
+        ix = z1 > 0;
+        ymax = max(sum(z1.*ix));
+        ix = z1 < 0;
+        ymin = min(sum(z1.*ix));
+        if ymax-ymin < 1e-6
+            continue
+        end
+        figure('Name',M_.endo_names(i_var(j),:));
+        ax=axes('Position',[0.1 0.1 0.6 0.8]);
+        axis(ax,[xmin xmax ymin ymax]);
+        plot(ax,x(2:end),z1(end,:),'k-','LineWidth',2)
+        hold on;
+        for i=1:gend
+            i_1 = i-1;
+            yp = 0;
+            ym = 0;
+            for k = 1:exo_nbr+1 
+                zz = z1(k,i);
+                if zz > 0
+                    fill([x(i) x(i) x(i+1) x(i+1)],[yp yp+zz yp+zz yp],k);
+                    yp = yp+zz;
+                else
+                    fill([x(i) x(i) x(i+1) x(i+1)],[ym ym+zz ym+zz ym],k);
+                    ym = ym+zz;
+                end
+                hold on;
+            end
+        end
+        plot(ax,x(2:end),z1(end,:),'k-','LineWidth',2)
+        hold off;
+
+        axes('Position',[0.75 0.1 0.2 0.8]);
+        axis([0 1 0 1]);
+        axis off;
+        hold on;
+        y1 = 0;
+        height = 1/(exo_nbr+1);
+        labels = strvcat(M_.exo_names,'Initial values');
+        
+        for j=1:exo_nbr+1
+            fill([0 0 0.2 0.2],[y1 y1+0.7*height y1+0.7*height y1],j);
+            hold on
+            text(0.3,y1+0.3*height,labels(j,:),'Interpreter','none');
+            hold on
+            y1 = y1 + height;
+        end
+        hold off
+    end

matlab/shock_decomposition.m

@@ -0,0 +1,73 @@
+function z = shock_decomposition(M_,oo_,varlist)
+% function z = shock_decomposition(R,epsilon,varlist)
+% Computes shocks contribution to a simulated trajectory
+%
+% INPUTS
+%    R:         mm*rr matrix of shock impact
+%    epsilon:   rr*nobs matrix of shocks
+%    varlist:   list of variables
+%
+% OUTPUTS
+%    z:         nvar*rr*nobs shock decomposition
+%
+% SPECIAL REQUIREMENTS
+%    none
+%  
+% part of DYNARE, copyright Dynare Team (2004-2008)
+% Gnu Public License.
+    
+% number of variables
+    endo_nbr = M_.endo_nbr;
+
+% number of shocks
+    nshocks = M_.exo_nbr;
+
+    % indices of endogenous variables
+    [i_var,nvar] = varlist_indices(varlist);
+
+    % reduced form
+    dr = oo_.dr;
+
+    % data reordering
+    order_var = dr.order_var;
+    inv_order_var = dr.inv_order_var;
+
+
+    % coefficients
+    A = dr.ghx;
+    B = dr.ghu;
+    
+    % initialization
+    for i=1:nshocks
+        epsilon(i,:) = eval(['oo_.SmoothedShocks.' M_.exo_names(i,:)]);
+    end
+    gend = size(epsilon,2);
+    
+    z = zeros(endo_nbr,nshocks+2,gend);
+    for i=1:endo_nbr
+        z(i,end,:) = eval(['oo_.SmoothedVariables.' M_.endo_names(i,:)]);
+    end
+
+    maximum_lag = M_.maximum_lag;
+    lead_lag_incidence = M_.lead_lag_incidence;
+    for i=1:gend
+        if i > 1 & i <= maximum_lag+1
+            lags = min(i-1,maximum_lag):-1:1;
+            k2 = dr.kstate(find(dr.kstate(:,2) <= min(i,maximum_lag)+1),[1 2]);
+            i_state = order_var(k2(:,1))+(min(i,maximum_lag)+1-k2(:,2))*M_.endo_nbr;
+        end
+        
+        if i > 1
+            tempx = permute(z(:,1:nshocks,lags),[1 3 2]);
+            m = min(i-1,maximum_lag);
+            tempx = [reshape(tempx,endo_nbr*m,nshocks); zeros(endo_nbr*(maximum_lag-i+1),nshocks)];
+            z(:,1:nshocks,i) = A(inv_order_var,:)*tempx(i_state,:);
+            lags = lags+1;
+        end
+
+        z(:,1:nshocks,i) = z(:,1:nshocks,i) + B(inv_order_var,:).*repmat(epsilon(:,i)',nvar,1);
+        z(:,nshocks+1,i) = z(:,nshocks+2,i) - sum(z(:,1:nshocks,i),2);
+    end
+    
+    
+