% This function states the model equations
% Output are analytical derivatives

% modified: KM 22/06/2005

function [fx,fxp,fy,fyp,fypyp,fypy,fypxp,fypx,fyyp,fyy,fyxp,fyx,fxpyp,fxpy,fxpxp,fxpx,fxyp,fxy,fxxp,fxx,f] = model

syms PAR_alpha PAR_beta PAR_delta PAR_psi PAR_gamma PAR_g PAR_mu PAR_Lbar PAR_taulbar PAR_taux PAR_rho
syms PAR_sigmatau PAR_sigmaz PAR_sigma
syms VAR_C VAR_K VAR_L VAR_Y VAR_PROD VAR_Cp VAR_Kp VAR_Lp VAR_Yp VAR_PRODp
syms VAR_X VAR_taul VAR_Xp VAR_taulp 


f1=VAR_C+(1+PAR_gamma)*VAR_Kp-(1-PAR_delta)/VAR_X*VAR_K-VAR_Y;
f2=(VAR_K/VAR_X)^(PAR_alpha)*VAR_L^(1-PAR_alpha)-VAR_Y;
f3=PAR_psi*VAR_L*((PAR_Lbar-VAR_L)^(-PAR_sigma))-(1-VAR_taul)*(1-PAR_alpha)*VAR_Y/VAR_C;
f4=PAR_beta*(((1-PAR_delta)*(1+PAR_taux)+PAR_alpha*VAR_Yp*VAR_Xp/VAR_Kp)/VAR_Xp)/VAR_Cp-(1+PAR_taux)/VAR_C;
f5=log(VAR_Xp)-PAR_mu;
f6=VAR_taulp-PAR_rho*VAR_taul-(1-PAR_rho)*PAR_taulbar;
f7=VAR_PROD-VAR_Y/VAR_L;

f = [f1;f2;f3;f4;f5;f6;f7];

x = [VAR_K VAR_X VAR_taul];
y = [VAR_Y VAR_C VAR_L VAR_PROD];
xp = [VAR_Kp VAR_Xp VAR_taulp];
yp = [VAR_Yp VAR_Cp VAR_Lp VAR_PRODp];

f = subs(f, [x,y,xp,yp], exp([x,y,xp,yp]));

[fx,fxp,fy,fyp,fypyp,fypy,fypxp,fypx,fyyp,fyy,fyxp,fyx,fxpyp,fxpy,fxpxp,fxpx,fxyp,fxy,fxxp,fxx]=anal_deriv(f,x,y,xp,yp);