/*
things to add:
-~~~something to get stationarity, like debt elastic interest rate

-~~fireign assets denominated in terms of tradables-check

- checkc the balance and istribution of profits of the intermediate goods sector
- fixed cost in market clearing
-profits and busget contraint(includes fixed cost also or just vairable cost)

// issues
-number of forward looking variables. right now the exoegnous process has been made forward looking

*/
@#define newcost = 1
//1 for the new cost structure where the interest rate factors enter the price indices and mc
@#define logvariables = 1
@#define model = 0
//0 for value added model, 1 for model with imported intermediates
@#define tf = 1
//0 for no trade finacne, 1 for model with trade finance
var omega lambda n_t n_n wage_t wage_n p_t p_n p_m r c k_t k_n u_t u_n  f_t f_n g_t g_n g_m g
p_g inv_t inv_n debt x m tot tb gdp    mc_j p_v  y_t vaexp r_v r_m imp_final


@#if tf
a_bar r_b  a_nd a_d n_d collat_min profit_min


@#endif

@#if logvariables
C TOT GDP X M
@#endif

;
varexo eps;
predetermined_variables k_t k_n debt;
parameters  rho_t rho_n  bbeta ssigma aalpha_t aalpha_n eta
wg_t wg_n wg_m    ddelta  a_t a_n  ubar_t ubar_n
p_tstar p_mstar rbar dbar psi totbar rho_tot mc_b,dwl_b,d_b pareto1 a_min mu
w_v w_m  a_av d_v d_m d_f v_f m_f
;



//CALIBRATING PARAMETERS
/*
w_v and w_m ar efor the production function in the tradable goods sector
p_v is the price of value added in th tradable sector
d_v and d_m: ractions that have to be financed with external borrowing
*/
@#if model
a_min = 0.103;
a_t =7.5;
@#else
a_min = 0.15;
//a_t =5;
a_t =4.2;
@#endif
a_n =1;

rho_t =1.45;
rho_n =1.45;
 

 
bbeta =0.97;
ssigma =2;
aalpha_t =0.33;
aalpha_n =0.33;



ddelta =0.025;
//a_t =2;



ubar_t =0.01;
ubar_n=0.01;
p_tstar=1;
 p_mstar=1;
rbar = (1/bbeta);

dbar =  0;
psi = 0.27;
rho_tot = 0.9;
totbar = 1;
dwl_b =1 ;
d_b =0.02 ;
eta = 3.8;
mc_b = (1/(1-dwl_b*d_b));
mu = eta/(eta-1);
pareto1 = 3.4;

@#if model
@#if tf

a_t = 12.8;
w_v = 0.1;
w_m = 0.5;
wg_t =0.2;
wg_n =0.7;
wg_m =0.1;
d_v = 0;
d_m = 0.5;
d_f = 0.5;
v_f = 1.45;
m_f = 0;
@#else
a_min = 0.7;
a_t = 1;
w_v = 0.65;
w_m = 1-w_v;
wg_t =0.4;
wg_n =0.52;
wg_m =1-wg_t-wg_n;
d_v = 0;
d_m = 0;
d_f = 0;
v_f = 1.4;
m_f = 0;
@#endif
@#else
@#if tf
//a_t = 1;
w_v = 1;
w_m = 0;
wg_t =0.3;
wg_n =0.6;
wg_m =1-wg_t-wg_n;
d_v = 0.5;
d_m = 0;
d_f = 0;
v_f = 1;
m_f = 0;
@#else
w_v = 1;
w_m = 0;
wg_t =0.3;
wg_n =0.6;
wg_m =1-wg_t-wg_n;
d_v = 0;
d_m = 0;
d_f = 0;
v_f = 1;
m_f = 0;
@#endif

@#endif

a_av = a_min*((pareto1/(pareto1-eta+1))^(1/(eta-1)));













rho = 0.9;






// nu is the parameter that goes with H in the utility funciton
// etac is the elasticity of substitution for consumption (b/w the three components)
// a_nt and a_t are prdoductivity numbers
/*
n_n, n_t = labor used in production of non tradable goods and tradable goods respectively 
x: exports
m: imports
xbar is the steady state value of exports, just calibrate it to some random value. it is not important as long as
you are not interested in the value of the debt variable
note that the variable debt is in levels not log deviation
*/
rho = 0.9;


model;

//First order conditions for household problem

omega = c-(((n_t)^(rho_t))/(rho_t))-(((n_n)^(rho_n))/(rho_n));
omega^(-ssigma) = lambda*p_g;

(omega^(-ssigma))*((n_t)^(rho_t-1))=lambda*wage_t;
(omega^(-ssigma))*((n_n)^(rho_n-1))=lambda*wage_n;

//Euler equation
lambda*p_t(+1) = bbeta*r*lambda(+1)*p_t;

//forst order conditins wrt investment
lambda*p_g*(1+k_t(+1)-k_t)/lambda(+1)=bbeta*(p_g(+1)*(1-ddelta+k_t(+2)-k_t(+1))   +u_t(+1));
lambda*p_g*(1+k_n(+1)-k_n)/lambda(+1)=bbeta*(p_g(+1)*(1-ddelta+k_n(+2)-k_n(+1))   +u_n(+1));

//Tradeble and non tradable goods producing firms
u_t*k_t = aalpha_t*p_v*a_t*f_t;
u_n*k_n = aalpha_n*p_n*a_n*f_n;

wage_t*n_t = (1-aalpha_t)*p_v*a_t*f_t;
wage_n*n_n = (1-aalpha_n)*p_n*a_n*f_n;

f_t = (k_t^aalpha_t)*(n_t^(1-aalpha_t));
f_n = (k_n^aalpha_n)*(n_n^(1-aalpha_n));


//Market clearing

@#if newcost
a_t*f_t=v_f+w_v*y_t*((a_av)^(eta-1))*((mu*p_v*r_v/p_t)^(-eta));
@#else
a_t*f_t=v_f+w_v*y_t*((a_av)^(eta-1))*((mu*p_v/p_t)^(-eta));
@#endif
a_n*f_n=g_n;
y_t = g_t+x;
//Zero profit conditions
//p_t*a_t*f_t=wage_t*n_t+u_t*k_t;
//p_n*a_n*f_n=wage_n*n_n+u_n*k_n;


//Final goods-home retailers
g_t = g*wg_t*((p_t/p_g)^(-eta));
g_n = g*wg_n*((p_n/p_g)^(-eta));
g_m = g*wg_m*((p_m/p_g)^(-eta));

@#if newcost
m = g_m+m_f+w_m*y_t*((a_av)^(eta-1))*((mu*r_m*p_m/p_t)^(-eta));
@#else
m = g_m+m_f+w_m*y_t*((a_av)^(eta-1))*((mu*p_m/p_t)^(-eta));
@#endif

//final goods price index definition
p_g = (wg_t*(p_t^(1-eta))+wg_n*(p_n^(1-eta))+wg_m*(p_m^(1-eta)))^(1/(1-eta));

//market clearing;
g = c+inv_t+inv_n+0.5*((k_t(+1)-k_t)^2)+0.5*((k_n(+1)-k_n)^2);

//Capital accumulation equations
k_t(+1) = (1-ddelta)*k_t+inv_t;
k_n(+1) = (1-ddelta)*k_n+inv_n;


//current account equation
p_t*x-p_m*m=p_t*((debt(+1)/r) -debt);
//debt = 0;

r = rbar+psi*(exp(debt-dbar)-1);

//gdp = c+inv_n+inv_t+p_t*x-p_m*m;
gdp = c+inv_t+inv_n+0.5*((k_t(+1)-k_t)^2)+0.5*((k_n(+1)-k_n)^2)+p_t*x-p_m*m;
//gdp = 1;
//p_t = p_tstar;
//p_m = p_mstar;


tb = p_t*x-p_m*m;


tot = p_t/p_m;

//exogenous pricess for tot
log(tot(+1)/totbar) = rho_tot*log(tot/totbar)+eps;
//tot = totbar;


//additional equations for financial accelerator model

//definitions: marginal cost component and interest rate factor (corresponding to M in watson)

@#if newcost
mc_j = (w_v*((r_v*p_v)^(1-eta))+w_m*((r_m*p_m)^(1-eta)))^(1/(1-eta));
@#else
mc_j = (w_v*(p_v^(1-eta))+w_m*(p_m^(1-eta)))^(1/(1-eta));
@#endif
p_t = mu*mc_j/a_av;

//definitions of average productivities:

@#if tf

a_nd = a_bar*((pareto1/(pareto1-eta+1))^(1/(eta-1)));
a_d =(a_av/a_min)*(((1/n_d)*(a_min^(eta-1))-((1-n_d)/n_d)*(a_bar^(eta-1)))^(1/(eta-1)));

n_d = 1-(a_min/a_bar)^pareto1;


//definitions
r_v=1-d_v+d_v*r_b;
r_m=1-d_m+d_m*r_b;





//bank participation contraint
/*
pc1+pc2 = pc3
*/
#pc11 = ((mu/p_t)^(-eta))*y_t*(a_nd^(eta-1));
#pc12 = d_v*w_v*(r_v^(-eta))*(p_v^(1-eta))+d_m*w_m*(r_m^(-eta))*(p_m^(1-eta));
#pc13 = d_v*p_v*v_f+d_m*p_m*m_f;
#pc1 = (1-n_d)*r_b*((pc11*pc12)+pc13);

#pc21 =  ((mu/p_t)^(-eta))*y_t*(a_d^(eta-1));
#pc22 = mu*((mc_j)^(1-eta))-w_v*(r_v^(-eta))*(p_v^(1-eta))-w_m*(r_m^(-eta))*(p_m^(1-eta));
#pc23 = p_v*v_f+p_m*m_f;

#pc2 = n_d*((pc21*pc22)-pc23);

#pc31 = ((mu/p_t)^(-eta))*y_t*(a_av^(eta-1));
#pc32 = d_v*w_v*(r_v^(-eta))*(p_v^(1-eta))+d_m*w_m*(r_m^(-eta))*(p_m^(1-eta));
#pc33 = d_v*p_v*v_f+d_m*p_m*m_f;

#pc3 = mc_b*(pc31*pc32+pc33);

pc1+pc2 = pc3;


//zero profit condition which determines cutoff productivity:
//(((mu^(-eta))*a_t*f_t)/(((a_bar)^(1-eta))*((p_t)^(-eta))))*(mu*(mc_j)^(1-eta)-w_v*r_v*(p_v^(1-eta))-w_m*r_m*(p_m^(1-eta))-v_f*p_v-m_f*p_m)=0;
@#if newcost
(((mu^(-eta))*y_t)/(((a_bar)^(1-eta))*((p_t)^(-eta))))*(mu*(mc_j)^(1-eta)-w_v*r_v*(r_v^(-eta))*(p_v^(1-eta))-w_m*r_m*(r_m^(-eta))*(p_m^(1-eta)))=r_v*v_f*p_v+r_m*m_f*p_m;
@#else
(((mu^(-eta))*y_t)/(((a_bar)^(1-eta))*((p_t)^(-eta))))*(mu*(mc_j)^(1-eta)-w_v*r_v*(r_v^(-eta))*(p_v^(1-eta))-w_m*r_m*(r_m^(-eta))*(p_m^(1-eta)))-v_f*p_v-m_f*p_m=0;

@#endif


//collateral of the minimum productivity firm, shoudl not be negative

@#if newcost
collat_min = (((mu^(-eta))*y_t)/(((a_min)^(1-eta))*((p_t)^(-eta))))*(mu*(mc_j)^(1-eta)-w_v*(r_v^(-eta))*(p_v^(1-eta))-w_m*(r_m^(-eta))*(p_m^(1-eta)))-v_f*p_v-m_f*p_m;

profit_min = (((mu^(-eta))*y_t)/(((a_min)^(1-eta))*((p_t)^(-eta))))*(mu*(mc_j)^(1-eta)-w_v*r_v*(r_v^(-eta))*(p_v^(1-eta))-w_m*r_m*(r_m^(-eta))*(p_m^(1-eta)))-v_f*p_v-m_f*p_m;

@#else
collat_min = (((mu^(-eta))*y_t)/(((a_min)^(1-eta))*((p_t)^(-eta))))*(mu*(mc_j)^(1-eta)-w_v*(p_v^(1-eta))-w_m*(p_m^(1-eta)))-v_f*p_v-m_f*p_m;

profit_min = (((mu^(-eta))*y_t)/(((a_min)^(1-eta))*((p_t)^(-eta))))*(mu*(mc_j)^(1-eta)-w_v*r_v*(p_v^(1-eta))-w_m*r_m*(p_m^(1-eta)))-v_f*p_v-m_f*p_m;

@#endif
@#else
r_v=1;
r_m = 1;
@#endif

p_g = 1;

vaexp = p_m*m/gdp;
imp_final = g_m/m;
//Definitions of log variables
@#if logvariables
C = log(c);
TOT = log(tot);
GDP = log(gdp);
X = log(x);
M = log(m);
@#endif




end;
initval;

tot = totbar;
p_t = 1;
p_m = p_t/totbar;
p_g = 1;

@#if tf
r_b = 1.1;

r_v=1-d_v+d_v*r_b;
r_m=1-d_m+d_m*r_b;
@#else
r_v = 1;
r_m = 1;
@#endif


u_t = (1/bbeta)+ddelta-1;
u_n = (1/bbeta)+ddelta-1;


p_n = ((p_g^(1-eta)-wg_t*p_t^(1-eta)-wg_m*p_m^(1-eta))/(wg_n))^(1/(1-eta));
//p_v = p_n;
p_v = (p_t/r_v)*((((a_av/mu)^(1-eta))-w_m*(r_m^(1-eta)))/w_v)^(1/(1-eta));

//(((a_av/mu)^(1-eta))-w_m*(r_m^(1-eta)))/w_v;


wage_t = (a_t*p_v*(aalpha_t^(aalpha_t))*((1-aalpha_t)^(1-aalpha_t))/(u_t^aalpha_t))^(1/(1-aalpha_t));

n_t = (wage_t/p_g)^(1/(rho_t-1));
k_t = (aalpha_t/(1-aalpha_t))*n_t*wage_t/u_t;
f_t = (k_t^(aalpha_t))*(n_t^(1-aalpha_t));
inv_t = ddelta*k_t;




wage_n = (a_n*p_n*(aalpha_n^(aalpha_n))*((1-aalpha_n)^(1-aalpha_n))/(u_n^aalpha_n))^(1/(1-aalpha_n));

n_n = (wage_n/p_g)^(1/(rho_n-1));
k_n = (aalpha_n/(1-aalpha_n))*n_n*wage_n/u_n;
f_n = (k_n^(aalpha_n))*(n_n^(1-aalpha_n));
inv_n = ddelta*k_n;

g_n = a_n*f_n;

g_t = g_n*(wg_t/wg_n)*(p_n/p_t)^(eta);
g_m = g_n*(wg_m/wg_n)*(p_n/p_m)^(eta);
g = ((((wg_t)^(1/eta))*(g_t)^((eta-1)/eta)) +(((wg_n)^(1/eta))*(g_n)^((eta-1)/eta)) 
+(((wg_m)^(1/eta))*(g_m)^((eta-1)/eta)))^(eta/(eta-1))  ;
c = g-inv_t-inv_n;
omega = c-(((n_t)^(rho_t))/(rho_t))-(((n_n)^(rho_n))/(rho_n));
lambda = (omega^(-ssigma))/(p_g);
r = rbar;
debt = dbar;





//p_v = p_n;





 


m = g_m*(1+(wg_t/wg_m)*((p_m/p_t)^(eta))*w_m*((a_av)^(eta-1))*((mu*r_m)^(-eta)))/(1-w_m*((a_av)^(eta-1))*((mu*r_m)^(-eta)));
x=m/tot+((dbar/r) -dbar);
y_t = x+g_t;
mc_j = (w_v*((r_v*p_v)^(1-eta))+w_m*((r_m*p_m)^(1-eta)))^(1/(1-eta));

@#if tf
a_bar = ((((mu^(-eta))*y_t)/(((1)^(1-eta))*((p_t)^(-eta))))*((mu*(mc_j)^(1-eta)-w_v*((r_v*p_v)^(1-eta))-w_m*((r_m*p_m)^(1-eta)))/(r_v*v_f*p_v+r_m*m_f*p_m)))^(1/(1-eta));

n_d = 1-(a_min/a_bar)^pareto1;

a_nd = a_bar*((pareto1/(pareto1-eta+1))^(1/(eta-1)));
a_d =(a_av/a_min)*(((1/n_d)*(a_min^(eta-1))-((1-n_d)/n_d)*(a_bar^(eta-1)))^(1/(eta-1)));

//m = g_m+m_f+w_m*y_t*((a_av)^(eta-1))*((mu*r_m)^(-eta));
//brac = ((a_av)^(eta-1))*((mu*r_m)^(-eta))


@#endif



gdp = c+inv_t+inv_n+p_t*x-p_m*m;

@#if tf
collat_min = (((mu^(-eta))*y_t)/(((a_min)^(1-eta))*((p_t)^(-eta))))*(mu*(mc_j)^(1-eta)-w_v*(r_v^(-eta))*(p_v^(1-eta))-w_m*(r_m^(-eta))*(p_m^(1-eta)))-v_f*p_v-m_f*p_m;

profit_min = (((mu^(-eta))*y_t)/(((a_min)^(1-eta))*((p_t)^(-eta))))*(mu*(mc_j)^(1-eta)-w_v*r_v*(r_v^(-eta))*(p_v^(1-eta))-w_m*r_m*(r_m^(-eta))*(p_m^(1-eta)))-v_f*p_v-m_f*p_m;
@#endif
//gdp = 1;
vaexp = p_m*m/gdp;
imp_final = g_m/m;


tb = 0;
@#if logvariables
C = log(c);
TOT = log(tot);
GDP = log(gdp);
X = log(x);
M = log(m);
@#endif
end;
resid;
steady(solve_algo=3,maxit=1000);
//write_latex_dynamic_model;

/*
shocks;
var eps =1;
end;

stoch_simul(irf=20,order=1) C TOT GDP p_t p_m X M;
*/