function usg_noise_model
%usg_noise_model computes analytically the 
%derivatives of the first-step equilibrium conditions of the  open economy model with stationary, nonstationary, and noise shocks and imperfect information presented in the book Open Economy Macroeconomics, by Martín Uribe and Stephanie Schmitt-Grohé (Princeton University Press) in  the chapter `Business Cycles in Emerging Countries: Productivity Shocks Versus Financial Frictions.'' 
%  The first-step equilibrium conditions take the form
%  E_t f(yp,y,xp,x) =0, 
%where x is a state vector and y a control vector, and 
%'p' denotes next-period.
%The derivatives of this 
%system evaluated at the deterministic steady state 
%are denoted fx, fxp, fy, and fyp Up to first order, the 
%dynamics of this system are governed by the system
%xp = hx * x 
%y = gx * x, 
%If y is observed with error, then the observable, o,  up 
%to first-order evolves according to
%o = y + w, 
%where Ew=0 and Ew*w' = varme, 
%and varme is a matrix of parameters. 
%
%Output: the file usg_noise_model_num_eval.m containing analytical expressions for f, fx, fxp, fy, fyp, and varme
%
%Calls: anal_deriv.m and anal_deriv_print2f.m
%
%© Martin Uribe and Stephanie Schmitt-Grohé, February 2, 2014

approx = 1; %Order of approximation desired

%Define parameters
syms BETTA SIGG GAMA DELTA ALFA PHI PSSI RHOZ RHOX RHON RSTAR DBAR G SIGMAZ SIGMAX	 SIGMAN 

%Define variables 
syms c cp k kp k1 k1p z zp zback zbackp  gxhat gxhatp  h hp d dp  yy yyp ivv ivvp tby tbyp cay cayp gy gc giv yyback cback ivvback gyp gcp givp yybackp cbackp ivvbackp ghat ghatp ghatback ghatbackp

syms eta1 eta2 eta3 eta4  eta6 
syms eta1p eta2p eta3p eta4p eta6p
syms s sp noise noisep
 
%Argument in expressions related to capital adjustment costs
ac = exp(G+ghat)*kp/k-exp(G);
acp = exp(G+ghatp)*k1p/k1-exp(G);

%Trade balance
tb = yy - c - ivv - PHI/2 * ac^2;

%Interest Rate
r = RSTAR + PSSI * (exp(dp-DBAR) - 1);

%Marginal utility of consumption
la = GAMA * c^(GAMA*(1-SIGG)-1) * (1-h)^((1-GAMA)*(1-SIGG));
lap = GAMA * cp^(GAMA*(1-SIGG)-1) * (1-hp)^((1-GAMA)*(1-SIGG));

%Write equations (e1, e2,...en)
e1 = d - tb - dp/(1+r)*exp(G+ghat);

e2 = -yy + exp((1-ALFA)*(G+ghat)) * k^ALFA *  h^(1-ALFA);

e3 = -ivv + kp*exp(G+ghat) - (1-DELTA) *k;

e4 = - la  + BETTA * (1+r) * exp((G+ghat)*(GAMA*(1-SIGG)-1)) * lap;

e5 = - (1-GAMA)/GAMA * c / (1-h)+ (1-ALFA) * exp((1-ALFA)*(G+ghat)) * (k/h)^ALFA;

e6 = -la * (1+ PHI * ac) + BETTA * exp((G+ghat)*(GAMA*(1-SIGG)-1))* lap * (1 - DELTA + ALFA * exp((1-ALFA)*(G+ghatp)) * (hp/kp)^(1-ALFA) + PHI * (acp +exp(G)) * acp - PHI/2 * acp^2 );

e7 = -k1 + kp;

%trade-balance-to-output ratio
e8 = -tby + tb / yy; 

%Output growth
e9 = -gy + yy/yyback*exp(G+ghatback);

%Consumpiton Growth
e10 = -gc + c/cback*exp(G+ghatback);

%Investment growth rate
e11 = -giv + ivv/ivvback*exp(G+ghatback);

e12 = -yybackp + yy;

e13 = -cbackp + c;

e14 = -ivvbackp + ivv;

e15 = -ghatbackp + ghat;

%Matrices of the Kalman filter
F = [
RHOZ 0 0 0
 1 0 0 0
0 0 RHOX 0
0 0 0 RHON];

Hp = [
1 -1 1 0
0 0 1 1];

syms k11 k12 k21 k22 k31 k32 k41 k42
K = [
k11 k12
k21 k22
k31 k32
k41 k42];

e16 = -[eta1p; eta2p; eta3p; eta4p] + (F-K*Hp)*[eta1; eta2; eta3; eta4] + K*[ghat;s];

e17 = -[ghatp;sp]  + Hp*[eta1p;eta2p;eta3p;eta4p];
%Create function f
f = [e1;e2;e3;e4;e5;e6;e7;e8;e9;e10;e11;e12;e13;e14;e15;e16;e17];

% Define the vector of controls, y, and states, x
statevar_cu = [   yyback cback ivvback ghatback k d  eta1 eta2 eta3 eta4  ghat s];
statevar_cup = [  yybackp cbackp ivvbackp ghatbackp kp dp   eta1p eta2p eta3p eta4p ghatp sp];

controlvar_cu = [gy gc giv tby h yy c ivv k1];

controlvar_cup = [gyp gcp givp tbyp hp yyp cp ivvp k1p];

log_linearize_cu = [yyback cback ivvback k d gy gc giv h yy c ivv k1];
log_linearize_cup = [yybackp cbackp ivvbackp kp dp gyp gcp givp hp yyp cp ivvp k1p];

%variables to substitute from levels to logs
log_linearize = [log_linearize_cu log_linearize_cup];

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

%Compute analytical derivatives of f
[fx,fxp,fy,fyp]=anal_deriv(f,statevar_cu,controlvar_cu,statevar_cup,controlvar_cup,approx);

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

%Eliminate future variables
cu = [statevar_cu controlvar_cu];
cup = [statevar_cup controlvar_cup];

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);

%what follows is a fake construction of ETASHOCK. Under noise shocks, the correct matrix is constructed at the end of the program gx_hx_inputs.m
syms SIGMAX
ns = length(statevar_cu);
ETASHOCK(ns,1) = SIGMAX;

varshock = ETASHOCK*ETASHOCK';

syms STDmey STDmec STDmeiv STDmetby
varme = diag([STDmey; STDmec; STDmeiv; STDmetby
].^2);

%Print derivatives to file `model_num_eval.m'  for model evaluation
filename = 'usg_noise_model';
anal_deriv_print2f(filename,fx,fxp,fy,fyp,f,ETASHOCK,varme);

