%mxns.m
%Produce the equilibrium conditions and their first-order derivatives  of the MXN model with interest-rate spreads developed in the paper ``How Important Are Terms Of Trade Shocks? by  S. Schmitt-Grohe and M. Uribe.
%Run this program only once or every time the structure of the model changes. No need to re-run it 
%when parameter values change. 
%This program produces a  log-linear approximation to the  function f, which defines  the DSGE model: 
%  E_t f(yp,y,xp,x) =0. 
%here, a p denotes next-period variables.  
%
%Output: Analytical expressions for f and its first derivatives as well as x and y. 
%The output is written to <filename>_num_eval.m which can then be run for numerical evaluations
%
%Calls: anal_deriv.m and first_order_print2f.m 
%
% ©  Stephanie Schmitt-Grohe and Martin Uribe, 2016. 
%
%Date September 2016
clear all 
cesm =0

filename = 'mxns'; %a file with this name and suffix _num_eval.m is created at the end of this program and contains symbolic expression for the function f and its derivatives fx fy fxp fyp 

syms OMEGAM OMEGAX SIGG PHIM PHIX  PHINT BETTA DELTA ALFAM ALFAX  RSTAR DBAR  PSSI TOT S OUTPUT TFPX TFPM TFPN
syms  STD_EPS_A ETASHOCK   RHO_TOTTOT RHO_TOTS RHO_STOT RHO_SS STD_TOTTOT STD_TOTS STD_STOT STD_SS
syms OMEGAN CHITAU CHIM MUTN MUMX ALFAN
syms PM PX PN

syms c cp hm hmp h_x h_xp d dp im imp ix ixp  km kmp kmfu kmfup kx kxp kxfu kxfup am amp ax axp px pxp pm pmp ym ymp yx yxp r rp wm wmp wx wxp ux uxp um ump  m mp x xp  la lap 
syms tot totp s sp tfp_x tfp_xp tfp_m tfp_mp  tfp_n tfp_np output outputp  tby tbyp ivv ivvp
syms hn hnp ivn ivnp kn knp wn wnp ptau ptaup un unp an anp yn ynp pn pnp  atau ataup ataun knfu knfup
syms output_constant_prices output_constant_pricesp c_constant_prices c_constant_pricesp ivv_constant_prices ivv_constant_pricesp

%Utitlity function: 
U =((c - hn^OMEGAN/OMEGAN - hm^OMEGAM/OMEGAM -  h_x^OMEGAX/OMEGAX)^(1-SIGG)-1)/(1-SIGG);


if cesm==1
A = (CHIM * am^(1-1/MUMX ) + (1-CHIM) * ax^(1-1/MUMX))^(1/(1-1/MUMX));
else
A = am^CHIM*ax^(1-CHIM);
end



B = (CHITAU * atau^(1-1/MUTN ) + (1-CHITAU) * an^(1-1/MUTN))^(1/(1-1/MUTN));





%Equilibrium conditions. The symbols e1, e2, ... denote equation 1, equation2, ...

%U_1 = la
e1 = -la + diff(U, 'c'); 

%U_2/U_1 = wm
e2 = wm + diff(U, 'hm')/diff(U,'c'); 

%U_3/U1 = wx
e3 =  [wx +  diff(U, 'h_x')/diff(U,'c');wn +  diff(U, 'hn')/diff(U,'c')];

%Euler d
e4 = -la*ptau + BETTA * (1+r) * lap*ptaup ; 

%Euler km
e5 = -la *(1+ PHIM*(kmp-km)) + BETTA*lap*  (ump + 1-DELTA +  PHIM * (kmfup-kmp));

%Euler km
e6 = [-la *(1+ PHIX*(kxp-kx)) + BETTA*lap * (uxp + 1-DELTA +  PHIX * (kxfup-kxp)); -la *(1+ PHINT*(knp-kn)) + BETTA*lap * (unp + 1-DELTA +  PHINT * (knfup-knp))];

%Evolution kmp
e7 =  -kmp + (1-DELTA)*km + im; 

%Evolution kxp and knp
e8 =  [-kxp + (1-DELTA) *kx +  ix; -knp + (1-DELTA) *kn +  ivn]; 

%Final good aggregation
e9 = [-pm + ptau* diff(A, 'am'); -px + ptau*diff(A, 'ax')]; 

%Final good aggregation
e10 = [-ptau + diff(B, 'atau'); -pn + diff(B,'an')];

%Production of Importables
e11 = - tfp_m* km^ALFAM*hm^(1-ALFAM)+ ym ; 

%Production of Importables
e12 = [- tfp_x* kx^ALFAX*h_x^(1-ALFAX)+ yx ; -tfp_n*kn^ALFAN*hn^(1-ALFAN)+ yn] ; 

%Foc Firms capital
e13 = [-um + pm *tfp_m *ALFAM*km^(ALFAM-1)*hm^(1-ALFAM); -un + pn* tfp_n*ALFAN*kn^(ALFAN-1)*hn^(1-ALFAN)]; 

%Foc Firms labor
e14 = -ux + px *tfp_x *ALFAX*kx^(ALFAX-1)*h_x^(1-ALFAX); 

%Foc Firms capital
e15 = [-wm + pm *tfp_m *(1-ALFAM)*km^ALFAM*hm^(-ALFAM); -wn + pn *tfp_n *(1-ALFAN)*kn^ALFAN*hn^(-ALFAN)]; 

%Foc Firms labor
e16 = -wx  +  px *tfp_x *(1-ALFAX)*kx^ALFAX*h_x^(-ALFAX); 

% Market clearing for final goods: 
e17 = - B + c + im +ix +ivn + PHIX/2*(kxp-kx)^2 + PHIM/2*(kmp-km)^2+PHINT/2*(knp-kn)^2; 

%Imports at final goods prices 
e18 = [pm*(-am+ym) + m; ; -A + atau; an - yn];

%Exports at final good prices
e19 = px*(-ax+yx) -x;

%Evolution of foreign debt
e20 = -ptau*d - m + ptau*dp /(1+r) + x; 

%Country premium
e21 = -r + RSTAR + s + PSSI * (exp(dp-DBAR) -1);

%Definition of Tot
e22 = -tot + px / pm; 

%Evolution of TOT, s, and TFP (TFP shocks are not explored in this paper)
e23 = [-log(totp/TOT) + RHO_TOTTOT * log(tot/TOT) +  RHO_TOTS * s;
-sp + RHO_STOT * log(tot/TOT) + RHO_SS * s;
    -log(tfp_xp/TFPX) + RHO_TOTTOT * log(tfp_x/TFPX);
    -log(tfp_mp/TFPM) + RHO_TOTTOT * log(tfp_m/TFPM);
    -log(tfp_np/TFPN) + RHO_TOTTOT * log(tfp_n/TFPN);
       ];

%make kfu=kp
e24 = [-kxfu+kxp; -kmfu + kmp; -knfu + knp];

%definition of output in terms of final goods and output at constant prices
e25 = [-output+px*yx+pm*ym+pn*yn; 
            -ivv+im+ix+ivn; 
            -output_constant_prices+PX*yx+PM*ym+PN*yn;
            -c_constant_prices+c/output*output_constant_prices; 
            -ivv_constant_prices+ivv/output*output_constant_prices];

%trade balance in final goods relative to steady state output
e26 = -tby + (x/output-m/output)*output_constant_prices/OUTPUT; 

%Create function f
f = eval(eval([e1;e2;e3;e4;e5;e6;e7;e8;e9;e10;e11;e12;e13;e14;e15;e16;e17;e18;e19;e20;e21;e22;e23;e24; e25; e26]));

f_nosub = f;

% Define the vector of controls in periods t and t+1, controlvar and controlvarp, and the vector of states in periods t and t+1, statevar and statevarp
%States to substitute from levels to logs
states_in_logs = [ km kx kn tfp_m tfp_x tfp_n tot];
states_in_logsp = [ kmp kxp knp tfp_mp tfp_xp tfp_np totp];
statevar = [d states_in_logs s];
statevarp = [dp  states_in_logsp sp];
 
controls_in_logs = [c      hm  h_x    hn   im   ix    ivn  am   ax an  atau       px  pm    pn   ptau  ym    yx  yn   wm   wx   wn   ux   um  un   m   x   output    la   kmfu   kxfu  knfu  ivv output_constant_prices c_constant_prices ivv_constant_prices];
controls_in_logsp = [cp  hmp h_xp hnp imp ixp  ivnp amp axp anp ataup   pxp pmp pnp ptaup ymp yxp ynp wmp wxp wnp uxp ump unp mp xp outputp  lap kmfup kxfup knfup   ivvp output_constant_pricesp c_constant_pricesp ivv_constant_pricesp];
controlvar = [controls_in_logs tby r]; 
controlvarp = [controls_in_logsp tbyp rp]; 

%Number of states
ns = length(statevar);

%Make f a function of the logarithm of the state and control vector

%variables to substitute from levels to logs
variables_in_logs = transpose([states_in_logs, controls_in_logs, states_in_logsp, controls_in_logsp]);

f = subs(f, variables_in_logs, exp(variables_in_logs));

approx = 1;

%Compute analytical derivatives of f
[fx,fxp,fy,fyp]=anal_deriv(f,statevar,controlvar,statevarp,controlvarp,approx);

%Make f and its derivatives a function of the level of its arguments rather than the log
f = subs(f, variables_in_logs, log(variables_in_logs));
fx = subs(fx, variables_in_logs, log(variables_in_logs));
fy = subs(fy, variables_in_logs, log(variables_in_logs));
fxp = subs(fxp, variables_in_logs, log(variables_in_logs));
fyp = subs(fyp, variables_in_logs, log(variables_in_logs));

%Symbolically evaluate f and its derivatives at the nonstochastic steady state (c=cp, etc.)
cu = transpose([statevar controlvar km kx kn ]);
cup = transpose([statevarp controlvarp kmfu kxfu knfu ]);

for subcb=1:2 %substitution must be run twice  in case the original system is a stochastic difference equation of order higher than one. For example, the current model features k_t, k_t+1, and k_t+2 and thus is a 2nd order difference equation

f = subs(f, cup,cu,0);
fx = subs(fx, cup,cu,0);
fy = subs(fy, cup,cu,0);
fxp = subs(fxp, cup,cu,0);
fyp = subs(fyp, cup,cu,0);
end



%Construct ETASHOCK matrix, which determines the var/cov of the forcing term of the system. Specifically, the state vector evolves over time according to 
%x_t+1 = hx x_t + ETASHOCK epsilon_t+1
PI_tot_s = [STD_TOTTOT STD_TOTS;
                   STD_STOT STD_SS];
ETASHOCK(ns-1:ns,1:2) = PI_tot_s;
ETASHOCK(1,1) = 0;

varshock = ETASHOCK*ETASHOCK';

%Print derivatives to file <filename>_num_eval.m'  for model evaluation
first_order_print2f(filename,fx,fxp,fy,fyp,f,ETASHOCK,statevar,controlvar); %this script can be downloaded from http://www.columbia.edu/~mu2166/1st_order/1st_order.htm

eval(['save ' filename '.mat'])