%SP.M
function [fx,fxp,fy,fyp,f] = sp
%This program computes a log-linear approximation to the  function f for an economy with catching up with the Joneses good by good. The model is described in ``Deep Habits''
%The function f defines  the DSGE model (a p denotes next-period variables) : 
%  E_t f(yp,y,xp,x) =0. 
%
%Inputs: theta_s theta_d theta_g (the degree of superficial, deep, and government habit)
%
%Output: Numerical first derivative of f
%
%Calls: anal_deriv.m num_eval.m dha_ss.m
%
%(c) Stephanie Schmitt-Grohe and Martin Uribe
%
%Date December 23, 2003

%Define parameters
syms SIGMA CHI GAMA BETTA DELTA ALFA ETA RHOZ RHOG RHOV GBAR PHI CSTAR IVSTAR GSTAR  

%Define variables 

syms c cp h hp w wp u up k kp iv ivp mu mup a ap g gp output outputp v vp 

%Give functional form for  production, and utility functions

%Utility function
util = ((c-CSTAR-log(v))^(1-SIGMA) -1) / (1-SIGMA) + GAMA * ((1-h)^(1-CHI)-1) / (1-CHI);
utilp = ((cp-CSTAR-log(vp))^(1-SIGMA) -1) / (1-SIGMA) + GAMA * ((1-hp)^(1-CHI)-1) / (1-CHI);

pf = a*k^(ALFA)*h^(1-ALFA);
pfp = ap * kp^(ALFA)*hp^(1-ALFA);

%Write equations (ei, i=1:n)

%Labor supply
e1 = w + diff(util,'h') / diff(util,'c');

%Euler equation
e2 = -diff(util,'c') + BETTA *  diff(utilp,'cp') * (1-DELTA+up);

%Resource constraint
e3 = -(pf - PHI) + c + iv + g;

%Evoluiton of capital
e4 = -kp + (1-DELTA) * k + iv - IVSTAR;

%Markup
e5 = -mu + diff(pf,'h') / w;

%nu
e6 = (1-ETA) * output  + ETA * (CSTAR+GSTAR+IVSTAR) + ETA / mu * (output-CSTAR-GSTAR-IVSTAR);


%MRT = w/u
e7 = -w / u + diff(pf,'h') / diff(pf,'k');

%Technology shock
e8 = log(ap) - RHOZ * log(a);

%Government Purchases shock
e9 = log(gp/GBAR) - RHOG * log(g/GBAR);

%Evolution of preference shock
e10 = - log(vp) + RHOV * log(v);

%Definition of output
e11 = -output + pf - PHI;


%Create function f
f = [e1;e2;e3;e4;e5;e6;e7;e8;e9;e10;e11];

% Define the vector of controls in period t, y, and t+1, yp, and the vector of states in period t, x, and t+1, xp

x = [k v g a];

xp = [kp vp gp ap];

y = [mu output w h c iv u ];

yp = [mup outputp wp hp cp ivp up];


%Make f a function of the logarithm of the state and control vector
f = subs(f, [x,y,xp,yp], exp([x,y,xp,yp]));


approx = 1; %Order of approximation desired
%Compute analytical derivatives of f
[fx,fxp,fy,fyp]=anal_deriv(f,x,y,xp,yp,approx);