% first block should contain 2 equation to define c and c_abr
% second blokc ok: three gov exp concepts (g g_h g_f - 3 equations)
% third block defines investment as a function of the price of capital d: should have 4 equations for (d, d_abr and i and i_abr)
% fourth block: allocation of demand: should contain 2x6=12 equation for the variables (c_h, c_f, i_h, i_f, y_h y_f and same for abr variables)
% production block: 2x5=10 equations k_bar k rk z y 
% price setting and marginal cost: pi_h, pi_f_abr, s, s_abr, p_h, p_f_abr, p p_abr and q (=9 equations - 9 var)
% policy rate (1eq = r)
 
% total 2+3+4+12+10+9=40 equations and variables (
% you can drop P_h P_f 
% missing equations for l and l_abr 
% what is tau ?
 
var 
c r pi g q c_abr y_h y_f c_h c_h_abr 
c_f c_f_abr i i_abr i_h i_h_abr i_f i_f_abr 
p_h p_f p g_h g_f k u k_bar k_abr u_abr 
k_bar_abr r_k r_k_abr l l_abr s s_abr pi_h pi_f pi_abr d d_abr y y_abr 
tau w w_abr;
 
varexo e;
 
parameters sigma rho_g n theta rho phi_h phi_f C_bar Y_bar I_bar K_bar 
G_bar L_bar phi_h_abr phi_f_abr eta sigma_a beta kappa zeta psi_nu nu
a delta psi_c psi_k tau_bar kappa_i mu
rho_r phi_y phi_pi phi_g alpha sigma_c;
 
sigma = 1;
eta = 2;
sigma_a = 0.01;
beta = 0.99;
rho = 0.8;
nu = 1;
n = 0.1;
a = 2/3;
theta = 7;
rho_g = 0.988;
alpha = 0.75;
delta = 0.025;
phi_h = 0.69;
phi_f = 1 - phi_h;
phi_h_abr = (n/(1-n))*phi_h;
phi_f_abr = 1 - phi_h_abr;
L_bar = 1/3; %arbitrary setting value of L
K_bar = 3.577;
Y_bar = L_bar^a*K_bar^(1-a); 
G_bar = 0.2;
C_bar = .2634;
I_bar = 1.554;
psi_nu = (a/nu)/((1-a)/nu + 1);
kappa = (1-alpha)*(1-alpha*beta)/alpha;
zeta = 1/(1+psi_nu*theta);
psi_c = a/((1-a)*(1/nu)+1);
psi_k = ((1+1/nu)*(1-a))/((1-a)*(1/nu)+1);
tau_bar = 0.27;
kappa_i = 2.5;
rho_r = 0.75;
phi_y = 0.85;
phi_pi = 1.5;
phi_g = 0.5;
mu = 1;
sigma_c = (-1/sigma*((1-a/mu*Y_bar/C_bar)*1/(1+1/nu))^(-1))^(-1);
 
 
model;
 
% this block should contain 2 equation to define c and c_abr
c = c(+1) - sigma*(r - pi(+1));
% I guess that q is defined elsewhere, so I think you have to choose between equation 2 or 3 but both together will be conflicting
c_abr = c_abr(+1) - sigma*(r - pi_abr(+1));
 
% ok three gov exp concepts - 3 equations
g = rho_g*g(-1) + e;
g_h = g*n;
g_f = g*(1-n);
 
% this block defines investment as a function of the price of capital d: should have 4 equations for (d, d_abr and i and i_abr)
% should the term in c be here ? 
d + kappa_i*(beta*(i(+1) - i) - (i - i(-1))) + 1/sigma_c*c = 0;     % yes because capital demand also depends on discounted present conspt. and in the paper the authors derive this equation in the appendix
d_abr + kappa_i*(beta*(i_abr(+1) - i_abr) - (i_abr - i_abr(+1))) + sigma_c*c_abr = 0;
% should this equation not contain the discount rate (or using the euler equation an expression in c and c(+1) instead of just an expression in c(+1))
d = beta*(1-delta)*d(+1) + (1-beta*(1-delta))*(r_k(+1) - 1/sigma_c*c(+1));
d_abr = beta*(1-delta)*d_abr(+1) + (1-beta*(1-delta))*(r_k_abr(+1) + pi(+1) - q - pi_abr(+1) - 1/sigma_c*c_abr(+1));
 
% allocation of demand: should contain 2x6=12 equation for the variables (c_h, c_f, i_h, i_f, y_h y_f and same for abr variables)
% i and i_arbr are determined by their foc, so here are two equation too much: given i and i_h => i_f follows
i = phi_h*i_h + phi_f*i_f;
i_abr = phi_h_abr*i_h_abr + phi_f_abr*i_f_abr;
i_h = i - eta*p_h;
i_h_abr = i_abr - eta*(p_h - q);
c_h = c - eta*p_h;
c_f = c - eta*p_f;

% why allocating total y if this is already done for the components of y (c, i, g) and y has been defined as the sum of the components
% also why would you need variable at time(+1) and at time(t) combined here?
y = c + i + g_h;
y_abr = c_abr + i_abr + g_f;
y = a*l + (1-a)*k;
y_abr = a*l_abr + (1-a)*k_abr;
y_h = phi_h*C_bar/Y_bar*c_h + phi_h*I_bar/Y_bar*i_h + phi_h_abr*(1-n)/n*C_bar/Y_bar*c_h_abr + phi_h_abr*(1-n)/n*I_bar/Y_bar*i_h_abr + g_h;
y_f = phi_f*C_bar/Y_bar*c_f + phi_f*I_bar/Y_bar*i_f + phi_f_abr*n/(1-n)*C_bar/Y_bar*c_f_abr + phi_f_abr*n/(1-n)*I_bar/Y_bar*i_f_abr + g_f;
 
% production block: 2x5=10 equations k_bar k rk z y 
% you are missing labor supply conditions
%non linear exprssions
k_bar = (1-delta)*k_bar(-1) + delta*i;
k_bar_abr = (1-delta)*k_bar_abr(-1) + delta*i_abr;
k = u + k_bar(-1);
k_abr = u_abr + k_bar_abr(-1);

sigma_a = r_k/u;
sigma_a = (r_k_abr-q)/u_abr;

l = y + (1-a)/a*r_k - (1-a)/a*s;        
l_abr = y_abr + (1-a)/a*r_k_abr - (1-a)/a*s_abr;     

r_k = a*l - a*k + s;
r_k_abr = a*l_abr - a*k_abr + s_abr;

w = (a-1)/a*l + (1-a)/a*k + s;
w_abr = (a-1)/a*l_abr + (1-a)/a*k_abr + s_abr;
 
% price setting and marginal cost: pi_h, pi_f_abr, s, s_abr, p_h, p_f_abr, p p_abr and q (=9 equations - 9 var)
% if s is the marginal cost, why not substituting this in here ?
pi_h = beta*pi_h(+1) + kappa*zeta*(psi_nu*y_h - p_h + 1/sigma*psi_c*c + psi_c*tau_bar*(1-tau_bar)*tau + psi_k*r_k);
pi_f = beta*pi_f(+1) + kappa*zeta*(psi_nu*y_f - p_f + 1/sigma*psi_c*c_abr + psi_c*tau_bar/(1-tau_bar)*tau + psi_k*r_k_abr);
s = (1/nu + 1 - a)*(y + (1-a)/a*r_k - (1-a)/a*s) - (1-a)*(y - r_k + s) + 1/sigma*c + tau_bar/(1-tau_bar)*tau;
s_abr = (1/nu + 1 - a)*(y_abr + (1-a)/a*r_k_abr - (1-a)/a*s_abr) - (1-a)*(y_abr + (1-a)/a*r_k_abr + s_abr) + 1/sigma*c_abr + tau_bar/(1-tau_bar)*tau;

% do you need all this ? you have inflation Phillips curves from which you can derive the price level (p_h=p_h(-1)+pi_h) and these can be used to deine other aggregate price levels and q. The capital P variables are not necessary probably.
p_h = p_h(-1) + pi_h;
p_f = p_f(-1) + pi_f;
phi_h*p_h + phi_f*p_f = p;

%average inflation
pi = phi_h*pi_h + phi_f*pi_f;
pi_abr = phi_h_abr*pi_h + phi_f_abr*pi_f;

% definition of q
phi_h_abr*p_h - phi_f_abr*p_f = q;
 
% policy rate (1eq = r)
r = rho_r*r(-1) + (1-rho)*(phi_pi*(n*pi_h + (1-n)*pi_f) + phi_y*(n*y + (1 - n)*y_abr) + phi_g*(n*g_h + (1 - n)*g_f));

% labour income tax rate (distortionary taxation)
tau = n/(1-n)*(p_h + g_h) + (1-n)/n*(p_f + g_f) - l - w;
  
end;
 
initval;
 
c = 0;
r = 0;
pi = 0;
g = 0;
q = 0;
c_abr = 0;
y_h = 0;
y_f = 0;
c_h = 0;
c_h_abr = 0;
c_f = 0;
c_f_abr = 0;
i = 0;
i_abr = 0;
i_h = 0;
i_h_abr = 0;
i_f = 0;
i_f_abr = 0;
p_h = 0;
p_f = 0;
p = 0;
g_h = 0;
g_f = 0;
k = 0;
u = 1;
k_bar = 0;
k_abr = 0;
u_abr = 1;
k_bar_abr = 0;
r_k = 1;
r_k_abr = 1;
l = 0;
l_abr = 0;
s = 0;
s_abr = 0;
pi_h = 0;
pi_f = 0;
pi_abr = 0;
d = 0;
d_abr = 0;
y = 0;
y_abr = 0;
tau= 0;
 
end;
 
steady(solve_algo=3,maxit=1000000);
 
check(qz_zero_threshold=1e-20);
 
resid(1);
 
shocks;
var e;
stderr 0.01;
end;
 
stoch_simul(order=1,IRF=25) c, c_h, c_f, y, y_abr, pi, r, p_h, p_f, i_h, i_f, k, k_abr;