function output=DecR(state,shocks,DR)
%%%%
% conditional on a state matrix this gets the output for the endogenous
% variables.

% shocks=s(2:end)-varphi*s(1:end-1);
% state=[s(1:end-1),zeros(size(s(1:end-1))),shocks];

%% DR is oo_.dr

N=length(DR.g_1);
T=size(shocks,1);
next=zeros(T,N);
nstate=size(DR.g_1,2);
next=zeros(T,N);

perm3=[];
for i=1:nstate
    for j=1:nstate
        for k=1:nstate
            if k>=j & j>=i
                perm3=[perm3,[i,j,k]'];
            end
        end
    end
end

perm2=[];
for i=1:nstate
    for j=1:nstate
            if j>=i
                perm2=[perm2,[i,j]'];
            end
    end
end

i=0;
NS=size(state,2);
for vec=perm3
    i=i+1;
    %% mult deals with the fact that the cross-derivatives are only
    %% stored once in the decision rule
    mult=3*sum(vec(2:end)~=vec(1:end-1)); mult=max(mult,1);
    if vec(1)<=NS; s1=state(vec(1)); else s1=shocks(:,vec(1)-NS); end
    if vec(2)<=NS; s2=state(vec(2)); else s2=shocks(:,vec(2)-NS); end
    if vec(3)<=NS; s3=state(vec(3)); else s3=shocks(:,vec(3)-NS); end
    next=bsxfun(@plus,next,mult*bsxfun(@times,s1.*s2.*s3,DR.g_3(:,i)'));
end

%% this does the quadratic term in the decision rule
i=0;
for vec=perm2
    i=i+1;
    mult=2*(vec(2)~=vec(1));mult=max(mult,1);
    if vec(1)<=NS; s1=state(vec(1)); else s1=shocks(:,vec(1)-NS); end
    if vec(2)<=NS; s2=state(vec(2)); else s2=shocks(:,vec(2)-NS); end
    next=bsxfun(@plus,next,mult*bsxfun(@times,s1.*s2,DR.g_2(:,i)'));
end

%% linear part
for i=1:nstate
    if i<=NS; s1=state(i); else s1=shocks(:,i-NS); end
    next=bsxfun(@plus,next,bsxfun(@times,s1,DR.g_1(:,i)'));
end

next=bsxfun(@plus,next,DR.g_0');
next=bsxfun(@plus,next,DR.ys(DR.order_var)');

output=next;
% U_next=next(:,DR.inv_order_var(1));
% c_next=next(:,DR.inv_order_var(2));
% n_next=next(:,DR.inv_order_var(5));
% pi_next=next(:,DR.inv_order_var(3));
%