function system_piecewise_linear_jacob_33_revised(theta_in::Array{Float64,1},params_in::Dict{String}{Float64},e_grid::Array{Float64,1},g_grid::Array{Float64,1},z_grid::Array{Float64,1},d_grid::Array{Float64,1},Rhat_grid::Array{Float64,1},yhat_grid::Array{Float64,1},EG_nodes::Array{Float64,1},ER_nodes::Array{Float64,1},EZ_nodes::Array{Float64,1},ED_nodes::Array{Float64,1},weight::Array{Float64,1})

#global beta eta tau gamma nu chi_h psi1 psi2 rho_r phi r pibar pi_ss y_ss gstar rho_g R_ss rho_z rho_d
#global EG_nodes ER_nodes EZ_nodes ED_nodes weight nQ
#global ncoefs

# Map coefficients
phi    = params_in["phi"];
tau    = params_in["tau"];
eta    = params_in["eta"];
psi1   = params_in["psi1"];
psi2   = params_in["psi2"];
beta   = params_in["beta"];
gamma  = params_in["gamma"];
pi_ss  = params_in["pi_ss"];
pibar  = params_in["pibar"];
gstar  = params_in["gstar"];
rho_z  = params_in["rho_z"];
rho_g  = params_in["rho_g"];
rho_r  = params_in["rho_r"];
rho_d  = params_in["rho_d"];
chi_h  = params_in["chi_h"];
nu     = params_in["nu"];
r      = params_in["r"];
y_ss   = params_in["y_ss"];
R_ss   = params_in["R_ss"];
pi_ss  = params_in["pi_ss"];
c_ss   = params_in["c_ss"];


nQ=length(EG_nodes)
ncoefs = 7;

theta_pi_1::Array{Float64,1}, theta_yy_1::Array{Float64,1},theta_cc_1::Array{Float64,1}, theta_pi_2::Array{Float64,1}, theta_yy_2::Array{Float64,1},theta_cc_2::Array{Float64,1},dd0::Float64,dd1::Float64,dd2::Float64,dd3::Float64,dd4::Float64,dd5::Float64 = get_coefs_plc(theta_in,params_in);
#theta_pi_1, theta_yy_1,theta_cc_1, theta_pi_2, theta_yy_2,theta_cc_2,dd0,dd1,dd2,dd3,dd4,dd5 = get_coefs_plc(theta_in,params_in);

# Number of unknowns
nTheta = ncoefs*3+3;

# Derivative of ZLB and non-ZLB coefficients to unknowns:
t::Array{Float64,1},hh::Array{Float64,2} = eval_hh_33_revised(theta_in,params_in);


mm  = zeros(ncoefs*3,nTheta);
mm[1:7,1:7]     = Matrix{Float64}(I,ncoefs,ncoefs);
mm[8:14,8:14]   = Matrix{Float64}(I,ncoefs,ncoefs);
mm[15:21,15:21] = Matrix{Float64}(I,ncoefs,ncoefs);


#-----------------------------------
# Housekeeping
#-----------------------------------
# Allocate memory for residuals
nGrid = length(e_grid);

res1 = zeros(nGrid);
res2 = zeros(nGrid);
res3 = zeros(nGrid);

dRR1_dj= zeros(nGrid,nTheta);
dRR2_dj= zeros(nGrid,nTheta);
dRR3_dj= zeros(nGrid,nTheta);

# Allocate memory for integrands
integrand_1 = zeros(nQ,1);
integrand_2 = zeros(nQ,1);

integrand_1_jacob = zeros(nQ,nTheta);
integrand_2_jacob = zeros(nQ,nTheta);

# weights for jacobian integration
big_weight = repeat(weight,1,nTheta);

#Other allocations
#dpihat_prime_dcoefs=zeros(24);
#dyhat_prime_dcoefs=zeros(24);
#dchat_prime_dcoefs=zeros(24);

dpihat_prime_dcoefs::Adjoint=zeros(24)';
dyhat_prime_dcoefs::Adjoint=zeros(24)';
dchat_prime_dcoefs::Adjoint=zeros(24)';

#dpihat_prime_dcoefs=zeros(24)';
#dyhat_prime_dcoefs=zeros(24)';
#dchat_prime_dcoefs=zeros(24)';


#println(typeof(dpihat_prime_dcoefs))

psi_use=zeros(7)
state_i=zeros(7)'
#state_i::RowVector=zeros(7)'
psi_prime=ones(7)
state_i_prime=ones(7)

printed=0

count_cons=zeros(nGrid)
count_cons_prime=zeros(nGrid)

## MAIN LOOP
for i=1:nGrid

    # Collocation point
    er   = e_grid[i];
    g    = g_grid[i];
    z    = z_grid[i];
    d    = d_grid[i];
    Rhat0= Rhat_grid[i];
    yhat0= yhat_grid[i];

    psi_use=[1.0 ,er,g,z,d,yhat0,Rhat0];
    state_i=[1.0, er, g, z, d, yhat0, Rhat0]';
    #println(typeof(state_i))
    Rbar = dd0 + dd1*er + dd2*g + dd3*z + dd4*d + dd5*yhat0;

    # Check ZLB
    if Rhat0>Rbar
        # Approx. Policies at current point
        yhat  = dot(psi_use,theta_yy_1);
        pihat = dot(psi_use,theta_pi_1);
        chat  = yhat-g

        Rhat   = (1-rho_r)*psi1*pihat + (1-rho_r)*psi2*(yhat-yhat0 + z) + rho_r*Rhat0 + er;

        # Transform to levels;
        ppi    = pi_ss*exp.(pihat);
        y      = y_ss*exp.(yhat);
        R      = R_ss*exp.(Rhat);
        c      = c_ss*exp.(chat);

        #-----------------
        # For Jacobian
        #-----------------
        # Derivative of decision rules  wrt coefficients

        dpihat_dcoefs = state_i*mm[1:7,:];
        dyhat_dcoefs  = state_i*mm[8:14,:];
        #dchat_dcoefs  = state_i*mm[15:21,:];
        dchat_dcoefs  = dyhat_dcoefs;

        #=
        dpihat_dcoefs = zeros(24)'
        dyhat_dcoefs = zeros(24)'
        dchat_dcoefs = zeros(24)'

        dpihat_dcoefs[1:7]  = state_i
        dyhat_dcoefs[8:14]  = state_i
        dchat_dcoefs[15:21] = state_i
        =#

        #println(size(dpihat_dcoefs))
        #println(size(dyhat_dcoefs))
        dRhat_dcoefs  = (1-rho_r)*psi1*dpihat_dcoefs + (1-rho_r)*psi2*dyhat_dcoefs;

        # Derivative of monetary policy rule wrt pihat and yhat
        dRhat_dpihat  = (1-rho_r)*psi1;
        dRhat_dyhat   = (1-rho_r)*psi2;

    else
        count_cons[i]=1
        # Approx. Policies at current point
        yhat  = dot(psi_use,theta_yy_2);
        pihat = dot(psi_use,theta_pi_2);
        chat  = yhat-g
        Rhat  = log(1/R_ss);

        # Transform to levels
        ppi   = pi_ss*exp.(pihat);
        y     = y_ss*exp.(yhat);
        R     = R_ss*exp.(Rhat);
        c     = c_ss*exp.(chat);

        #-----------------
        # For Jacobian
        #-----------------

        # Derivative of decision rules  wrt coefficients
        dpihat_dcoefs = state_i*hh[1:7,:];
        dyhat_dcoefs  = state_i*hh[8:14,:];
        #dchat_dcoefs  = state_i*hh[15:21,:];
        dchat_dcoefs  = dyhat_dcoefs;
        dRhat_dcoefs  = zeros(1,24);

        # Derivative of monetary policy rule wrt pihat and yhat
        dRhat_dpihat  = 0.0;
        dRhat_dyhat   = 0.0;

    end

    # Compute endogenous objects
    ff = c^-tau*y*( (1-nu) - chi_h*(c^tau)*(y^(1/eta)) - 0.5*phi*(ppi-pibar)^2 + nu*phi*(ppi-pibar)*ppi );

    # For Jacobian: derivative of X' wrt coefficients
    dXp_dcoefs      = zeros(7,24);
    dXp_dcoefs[6,:] = dyhat_dcoefs;
    dXp_dcoefs[7,:] = dRhat_dcoefs;

    # Derivatives of endogenous objects wrt decision rules
    dC_dyhat        = 0.0;
    dC_dpihat       = 0.0;
    dC_dchat        = c_ss*exp(chat);

    dC_dyhat_prime  = 0.0;
    dC_dpihat_prime = 0.0;
    dC_dchat_prime  = 0.0;

    dff_dpihat =  (c^-tau)*y*(-phi*(ppi-pibar) + 2*nu*phi*ppi - nu*phi*pibar)*pi_ss*exp(pihat);
    dff_dyhat  =  y_ss*exp(yhat)*c^-tau*(1 - nu - 0.5*phi*(ppi - pibar)^2 + phi*nu*(ppi - pibar)*ppi) - (1+1/eta)*(y^(1/eta))*y_ss*exp(yhat);
    dff_dchat  =  -tau*c_ss*exp(chat)*c^(-tau-1)*y*(1 - nu - 0.5*phi*(ppi-pibar)^2 + nu*phi*(ppi-pibar)*ppi );

    dff_dyhat_prime   = 0.0;
    dff_dpihat_prime  = 0.0;
    dff_dchat_prime   = 0.0;

    psi_prime[6]=yhat
    psi_prime[7]=Rhat

    state_i_prime[6]=yhat
    state_i_prime[7]=Rhat

    # Integrals
    for q=1:nQ

        # Get g_prime at current node
        g1::Float64  =  rho_g*g + EG_nodes[q];
        z1::Float64  =  rho_z*z + EZ_nodes[q];
        d1::Float64  =  rho_d*d + ED_nodes[q];
        er1::Float64 =  ER_nodes[q];

        #psi_prime[1:5]=[1.0, er1::Float64, g1::Float64, z1::Float64, d1::Float64];
        psi_prime[2]=er1
        psi_prime[3]=g1
        psi_prime[4]=z1
        psi_prime[5]=d1

        state_i_prime[2]=er1
        state_i_prime[3]=g1
        state_i_prime[4]=z1
        state_i_prime[5]=d1

        Rbar_prime = dd0 + dd1*er1 + dd2*g1 +dd3*z1 + dd4*d1 + dd5*yhat;

        # Check ZLB
        if Rhat > Rbar_prime

            # Approximate t+1 policies
            yhat_prime  = dot(psi_prime,theta_yy_1);   # No ZLB
            pihat_prime = dot(psi_prime,theta_pi_1);   # No ZLB
            #chat_prime  = dot(psi_prime,theta_cc_1);
            chat_prime  = yhat_prime-g1

            # Transform to levels
            ppi_prime    = pi_ss*exp.(pihat_prime);
            y_prime      = y_ss*exp.(yhat_prime);
            c_prime      = c_ss*exp.(chat_prime);


            #--------------------------------------------------------------
            # These are the objects that should be computed with the
            # symbolic toolbox
            #-------------------------------------------------------------



            #=
            for i_t=1:24
                dpihat_prime_dcoefs[i_t] = theta_pi_1[1]*dXp_dcoefs[1,i_t]
                dyhat_prime_dcoefs[i_t] = theta_yy_1[1]*dXp_dcoefs[1,i_t]
                dchat_prime_dcoefs[i_t] = theta_cc_1[1]*dXp_dcoefs[1,i_t]
                for j_t=2:7
                dpihat_prime_dcoefs[i_t] =+ theta_pi_1[j_t]*dXp_dcoefs[j_t,i_t]
                dyhat_prime_dcoefs[i_t]  =+ theta_yy_1[j_t]*dXp_dcoefs[j_t,i_t]
                dchat_prime_dcoefs[i_t]  =+ theta_cc_1[j_t]*dXp_dcoefs[j_t,i_t]
                end
            end
            =#

            dpihat_prime_dcoefs = theta_pi_1'*dXp_dcoefs;
            dyhat_prime_dcoefs  = theta_yy_1'*dXp_dcoefs;
            #dchat_prime_dcoefs  = theta_cc_1'*dXp_dcoefs;


            for ii=1:ncoefs
                @inbounds dpihat_prime_dcoefs[ii]    += state_i_prime[ii]
                @inbounds dyhat_prime_dcoefs[ii+7]   += state_i_prime[ii]
                #@inbounds dchat_prime_dcoefs[ii+14] += state_i_prime[ii]
            end

            dchat_prime_dcoefs  = dyhat_prime_dcoefs;

            else
            count_cons_prime[i]=count_cons_prime[i]+1
            # Approximate t+1 policies
            yhat_prime   = dot(psi_prime,theta_yy_2);
            pihat_prime  = dot(psi_prime,theta_pi_2);
            #chat_prime   = dot(psi_prime,theta_cc_2);
            chat_prime  = yhat_prime-g1

            # Transform to levels
            ppi_prime    = pi_ss*exp.(pihat_prime);
            y_prime      = y_ss*exp.(yhat_prime);
            c_prime      = c_ss*exp.(chat_prime);

            #  R_prime     = log(1/R_ss);

            #--------------------------------------------------------------
            # These are the objects that should be computed with the
            # symbolic toolbox
            #-------------------------------------------------------------

            #=
            dpihat_prime_dcoefs = theta_pi_2'*dXp_dcoefs;
            dyhat_prime_dcoefs  = theta_yy_2'*dXp_dcoefs;
            dchat_prime_dcoefs  = theta_cc_2'*dXp_dcoefs;


            for ii=1:24
                @inbounds dpihat_prime_dcoefs[ii]   = dot(state_i_prime,hh[1:7,ii])+dot(theta_pi_2,dXp_dcoefs[:,ii])
                @inbounds dyhat_prime_dcoefs[ii]    = dot(state_i_prime,hh[8:14,ii])+dot(theta_yy_2,dXp_dcoefs[:,ii])
                @inbounds dchat_prime_dcoefs[ii]    = dot(state_i_prime,hh[15:21,ii])+dot(theta_cc_2,dXp_dcoefs[:,ii])
            end
            =#

            dpihat_prime_dcoefs = state_i_prime'*hh[1:7,:]   + theta_pi_2'*dXp_dcoefs;
            dyhat_prime_dcoefs  = state_i_prime'*hh[8:14,:]  + theta_yy_2'*dXp_dcoefs;
            #dchat_prime_dcoefs = state_i_prime'*hh[15:21,:] + theta_cc_2'*dXp_dcoefs;
            dchat_prime_dcoefs  = dyhat_prime_dcoefs;

            #=
            dpihat_prime_dcoefs = theta_pi_2'*dXp_dcoefs;
            dyhat_prime_dcoefs  = theta_yy_2'*dXp_dcoefs;
            dchat_prime_dcoefs  = theta_cc_2'*dXp_dcoefs;

            for ii=1:ncoefs
                dpihat_prime_dcoefs[ii]    += state_i_prime[ii]
                dyhat_prime_dcoefs[ii+7]   += state_i_prime[ii]
                dpihat_prime_dcoefs[ii+14] += state_i_prime[ii]
            end

            if printed==0
            println(typeof(dpihat_prime_dcoefs))
            printed=1
            end
            =#
        end


        # Compute integrands at current node
        @inbounds integrand_1[q] =  beta*R.*(exp(d1-d)./(gamma.*exp(z1))).*(c_prime.^-tau).*(1 ./ppi_prime);
        @inbounds integrand_2[q] =  phi*nu*beta*(c_prime.^-tau).*exp(d1-d).*(y_prime).*(ppi_prime - pibar).*ppi_prime;

        #------------------------------
        # For Jacobian integration
        #------------------------------

        dY1_dpihat      = 0.0;
        dY1_dpihat_prime= 0.0;
        dY1_dyhat       = 0.0;
        dY1_dyhat_prime = y_ss*exp(yhat_prime);
        dY1_dchat       = 0.0;
        dY1_dchat_prime = 0.0;

        dR_dyhat         = R_ss*exp(Rhat)*dRhat_dyhat;
        dR_dpihat        = R_ss*exp(Rhat)*dRhat_dpihat;
        dR_dyhat_prime   = 0.0;
        dR_dpihat_prime  = 0.0;
        dR_dchat         = 0.0;
        dR_dchat_prime   = 0.0;

        dP1_dyhat        = 0.0;
        dP1_dpihat       = 0.0;
        dP1_dyhat_prime  = 0.0;
        dP1_dpihat_prime = pi_ss*exp(pihat_prime);
        dP1_dchat        = 0.0;
        dP1_dchat_prime  = 0.0;

        dC1_dpihat       = 0.0;
        dC1_dpihat_prime = 0.0;
        dC1_dyhat        = 0.0;
        dC1_dyhat_prime  = 0.0;
        dC1_dchat        = 0.0;
        dC1_dchat_prime  = c_ss*exp(chat_prime);

        # Derivative of Residuals for Jacobian
        dF1_dpihat       = -tau*c^(-tau-1)*dC_dpihat       - integrand_1[q]*((1/R)*dR_dpihat -tau/c_prime*dC1_dpihat -(1/ppi_prime)*dP1_dpihat);

        dF1_dpihat_prime = -tau*c^(-tau-1)*dC_dpihat_prime - integrand_1[q]*((1/R)*dR_dpihat_prime - tau/c_prime*dC1_dpihat_prime - (1/ppi_prime)*dP1_dpihat_prime);
        dF1_dyhat        = -tau*c^(-tau-1)*dC_dyhat        - integrand_1[q]*((1/R)*dR_dyhat - tau/c_prime*dC1_dyhat - (1/ppi_prime)*dP1_dyhat);
        dF1_dyhat_prime  = -tau*c^(-tau-1)*dC_dyhat_prime  - integrand_1[q]*((1/R)*dR_dyhat_prime  - tau/c_prime*dC1_dyhat_prime  - (1/ppi_prime)*dP1_dyhat_prime);
        dF1_dchat        = -tau*c^(-tau-1)*dC_dchat        - integrand_1[q]*((1/R)*dR_dchat  - tau/c_prime*dC1_dchat  - (1/ppi_prime)*dP1_dchat);
        dF1_dchat_prime  = -tau*c^(-tau-1)*dC_dchat_prime  - integrand_1[q]*((1/R)*dR_dchat_prime  - tau/c_prime*dC1_dchat_prime  - (1/ppi_prime)*dP1_dchat_prime);

        dF2_dpihat       = dff_dpihat       - integrand_2[q]*((-tau/c_prime)*dC1_dpihat + (1/y_prime)*dY1_dpihat + (1/(ppi_prime-pibar) + 1/ppi_prime)*dP1_dpihat);
        dF2_dpihat_prime = dff_dpihat_prime - integrand_2[q]*((-tau/c_prime)*dC1_dpihat_prime + (1/y_prime)*dY1_dpihat_prime + (1/(ppi_prime-pibar)+1/ppi_prime)*dP1_dpihat_prime);
        dF2_dyhat        = dff_dyhat        - integrand_2[q]*((-tau/c_prime)*dC1_dyhat + (1/y_prime)*dY1_dyhat + (1/(ppi_prime-pibar) + 1/ppi_prime)*dP1_dyhat);
        dF2_dyhat_prime  = dff_dyhat_prime  - integrand_2[q]*((-tau/c_prime)*dC1_dyhat_prime  + (1/y_prime)*dY1_dyhat_prime  + (1/(ppi_prime-pibar)+1/ppi_prime)*dP1_dyhat_prime);
        dF2_dchat        = dff_dchat        - integrand_2[q]*((-tau/c_prime)*dC1_dchat + (1/y_prime)*dY1_dchat + (1/(ppi_prime-pibar) + 1/ppi_prime)*dP1_dchat);
        dF2_dchat_prime  = dff_dchat_prime  - integrand_2[q]*((-tau/c_prime)*dC1_dchat_prime  + (1/y_prime)*dY1_dchat_prime  + (1/(ppi_prime-pibar)+1/ppi_prime)*dP1_dchat_prime);


        # THESE ARE THE INTEGRANDS FOR THE JACOBIAN

        for i_t=1:24
            @inbounds integrand_1_jacob[q,i_t] = dF1_dpihat*dpihat_dcoefs[i_t] + dF1_dpihat_prime*dpihat_prime_dcoefs[i_t] + dF1_dyhat*dyhat_dcoefs[i_t] + dF1_dyhat_prime*dyhat_prime_dcoefs[i_t] + dF1_dchat*dchat_dcoefs[i_t] + dF1_dchat_prime*dchat_prime_dcoefs[i_t];
            @inbounds integrand_2_jacob[q,i_t] = dF2_dpihat*dpihat_dcoefs[i_t] + dF2_dpihat_prime*dpihat_prime_dcoefs[i_t] + dF2_dyhat*dyhat_dcoefs[i_t] + dF2_dyhat_prime*dyhat_prime_dcoefs[i_t] + dF2_dchat*dchat_dcoefs[i_t] + dF2_dchat_prime*dchat_prime_dcoefs[i_t];
        end

    end

    dF3_dpihat       = phi*(ppi-pibar)*y*pi_ss*exp(pihat);
    dF3_dpihat_prime = 0.0;
    dF3_dyhat        = -(1/((gstar)*exp(g)) - 0.5*phi*(ppi-pibar)^2)*y_ss*exp(yhat);
    dF3_dyhat_prime  = 0.0;
    dF3_dchat        = c_ss*exp(chat);
    dF3_dchat_prime  = 0.0;


    # Computing Expectations
    #sol_int1 = integrand_1' * weight;
    #sol_int2 = integrand_2' * weight;

    sol_int1 = dot(integrand_1,weight);
    sol_int2 = dot(integrand_2,weight);


    # Computing Expectations for Jacobian

    #sol_int1_jacob = integrand_1_jacob'*weight;
    #sol_int2_jacob = integrand_2_jacob'*weight;

    sol_int1_jacob = sum(integrand_1_jacob.*big_weight,dims=1);
    sol_int2_jacob = sum(integrand_2_jacob.*big_weight,dims=1);

    # Storing Raw Jacobian
    dRR1_dj[i,:] = sol_int1_jacob;
    dRR2_dj[i,:] = sol_int2_jacob;
    #println(size(dF3_dpihat*dpihat_dcoefs))
    #println(size(dF3_dpihat_prime*dpihat_prime_dcoefs))
    #for i_t=1:24
    #    dRR3_dj[i,i_t] = dF3_dpihat*dpihat_dcoefs[i_t] + dF3_dpihat_prime*dpihat_prime_dcoefs[i_t] + dF3_dyhat*dyhat_dcoefs[i_t] + dF3_dyhat_prime*dyhat_prime_dcoefs[i_t] + dF3_dchat*dchat_dcoefs[i_t] + dF3_dchat_prime*dchat_prime_dcoefs[i_t];
    #end
    #dRR3_dj[i,:] = dF3_dpihat*dpihat_dcoefs + dF3_dpihat_prime*dpihat_prime_dcoefs + dF3_dyhat*dyhat_dcoefs + dF3_dyhat_prime*dyhat_prime_dcoefs + dF3_dchat*dchat_dcoefs + dF3_dchat_prime*dchat_prime_dcoefs;

    # Residual Function
    res1[i] = c^-tau - sol_int1[1];
    res2[i] = ff - sol_int2[1];
    #res3[i] = c - (1/((gstar)*exp(g)) - 0.5*phi*(ppi-pibar)^2)*y;

end

# Residuals to be minimized
#G = [res1;res2;res3];
G = [res1;res2];

#JACOB = [dRR1_dj; dRR2_dj; dRR3_dj];
JACOB= [dRR1_dj; dRR2_dj];
JACOB[:,15:21].=0
JACOB[:,24].=0

return G,JACOB,count_cons,count_cons_prime

end
