%   This is the .mod file to solve the deterministic model

%--------------------------{ Preamble }------------------------------------

%   This is where we establish the variables, exogenous variables, and the
%   parameters of the model.

var c c_p ca n n_p na m lam1 lam2 tau w pi d r r_d r_l h h_p ha y g l; 
varexo mu;
parameters xi theta pis delta phi alpha omega chi eta lambda taus beta x_d 
           rho_r rho_pi phi_g phi_n gs nas rs cst cpst cast nst npst nast 
           mst lam1st lam2st wst dst hst hpst hast yst gst lst cen cpen 
           caen nen npen naen men lam1en lam2en wen den hen hpen haen yen 
           gen len must muen;

load steady_state_values.mat

set_param_value('xi', xi);
set_param_value('theta', theta);
set_param_value('pis', pist);
set_param_value('delta', delta);
set_param_value('phi', phi);
set_param_value('alpha', alpha);
set_param_value('omega', omega);
set_param_value('chi', chi);
set_param_value('eta', eta);
set_param_value('lambda', lambda);
set_param_value('taus', taust);
set_param_value('beta', beta);
set_param_value('x_d', x_d);
set_param_value('rho_r', rho_r);
set_param_value('rho_pi', rho_pi);
set_param_value('phi_g', phi_g);
set_param_value('phi_n', phi_n);
set_param_value('gs', gst);
set_param_value('nas', nast);
set_param_value('rs', rst);
set_param_value('cst', cst);
set_param_value('cpst', cpst);
set_param_value('cast', cast);
set_param_value('nst', nst);
set_param_value('npst', npst);
set_param_value('nast', nast);
set_param_value('mst', mst);
set_param_value('lam1st', lam1st);
set_param_value('lam2st', lam2st);
set_param_value('wst', wst);
set_param_value('dst', dst);
set_param_value('hst', hst);
set_param_value('hpst', hpst);
set_param_value('hast', hast);
set_param_value('yst', yst);
set_param_value('gst', gst);
set_param_value('lst', lst);
set_param_value('must', must);
set_param_value('cen', cen);
set_param_value('cpen', cpen);
set_param_value('caen', caen);
set_param_value('nen', nen);
set_param_value('npen', npen);
set_param_value('naen', naen);
set_param_value('men', men);
set_param_value('lam1en', lam1en);
set_param_value('lam2en', lam2en);
set_param_value('wen', wen);
set_param_value('den', den);
set_param_value('hen', hen);
set_param_value('hpen', hpen);
set_param_value('haen', haen);
set_param_value('yen', yen);
set_param_value('gen', gen);
set_param_value('len', len);
set_param_value('muen', muen);



%----------------------------{ Model }-------------------------------------

%   Since this is a deterministic model, it does not have to be linearized.
%   Thus, we can simply write it in its nonlinear form.

model;
1/c - (c/m)^(chi-1)*(1/m) = lam2*tau;
n = xi*mu^(theta-1)*m*(lam1/(lam2-beta*lam2(+1)*(1-delta)/pi(+1)))^theta;
d = (1-xi)*m*(lam1/(lam2-beta*lam2(+1)*(r_d/pi(+1))))^theta;
1 = beta*(lam2(+1)/lam2)*(r/pi(+1));
1 = beta*(lam2(+1)/lam2)*(r_l/pi(+1));
eta = lam2*w;
(c/m)^chi = lam1*m;
tau*c_p = w*h_p;
w*h_p = mu*n_p;
h_p = 1/eta;
r_d = (1-omega)*r_l + (omega-x_d)*pi(+1)/(beta*lam2(+1)/lam2);
w = alpha*y/ha;
m = (xi^(1/theta)*(mu*n)^((theta-1)/theta)+(1-xi)^(1/theta)*d^((theta-1)/theta))^(theta/(theta-1));
tau*c + d + n + r_l(-1)/pi*l(-1) = w*h + r_d(-1)/pi*d(-1) + (1-delta)/pi*n(-1) + l;
l = (1-omega)*d;
y = ha^alpha*g^(1-alpha);
(tau-1)*ca + (1-phi)*(na-((1-delta)/pi)*na(-1)) = g;
log(tau/taus) = phi_g*log(g/gs) + phi_n*log((na-(1-delta)*na(-1))/(delta*nas));
ca = (1-lambda)*c+lambda*c_p;
ha = (1-lambda)*h+lambda*h_p;
na = (1-lambda)*n+lambda*n_p;
log(r/rs) = rho_r*log(r(-1)/rs)+rho_pi*log(pi(-1)/pis);
end;

%--------------------------{ Steady State Values }-------------------------

initval;
c=cst;
c_p=cpst; 
ca =cast;
n =nst;
n_p =npst;
na =nast;
m =mst;
lam1=lam1st; 
lam2 =lam2st;
tau =taust;
w =wst;
pi =pist;
d =dst;
r =rst;
r_d =(1-x_d)*rst;
r_l =rst;
h =hst;
h_p =hpst;
ha =hast;
y =yst;
g =gst;
l=lst;
mu=must;
end;
steady;

resid;

check;


endval;
c=cen;
c_p=cpen; 
ca =caen;
n =nen;
n_p =npen;
na =naen;
m =men;
lam1=lam1en; 
lam2 =lam2en;
tau =tauen;
w =wen;
pi =pien;
d =den;
r =ren;
r_d =(1-x_d)*ren;
r_l =ren;
h =hen;
h_p =hpen;
ha =haen;
y =yen;
g =gen;
l=len;
mu=muen;
end;
steady;


shocks;
var mu;
periods 1:8;
values 0;
end;

simul(periods=80);