% This program finds the policy rules of the model

% modified: KM 22/06/2005

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

[PAR_alpha, PAR_beta, PAR_delta, PAR_psi, PAR_gamma, PAR_g, PAR_mu,PAR_Lbar, PAR_taulbar, PAR_taux, PAR_rho, PAR_sigmatau, PAR_sigmaz, PAR_sigma]=parameters;

[VAR_C,VAR_K,VAR_L,VAR_Y,VAR_PROD, VAR_Cp,VAR_Kp,VAR_Lp,VAR_Yp, VAR_PRODp] = SSfinder;

VAR_X=PAR_mu;
VAR_Xp=PAR_mu;
VAR_taul=log(PAR_taulbar);
VAR_taulp=log(PAR_taulbar);

eta=[0 0; 1 0; 0 1];

%Order of approximation desired 
approx = 2;

%Obtain numerical derivatives of f
num_eval

%First-order approximation
[gx,hx] = gx_hx(nfy,nfx,nfyp,nfxp)

%Second-order approximation
[gxx,hxx] = gxx_hxx(nfx,nfxp,nfy,nfyp,nfypyp,nfypy,nfypxp,nfypx,nfyyp,nfyy,nfyxp,nfyx,nfxpyp,nfxpy,nfxpxp,nfxpx,nfxyp,nfxy,nfxxp,nfxx,hx,gx) 

[gss,hss] = gss_hss(nfx,nfxp,nfy,nfyp,nfypyp,nfypy,nfypxp,nfypx,nfyyp,nfyy,nfyxp,nfyx,nfxpyp,nfxpy,nfxpxp,nfxpx,nfxyp,nfxy,nfxxp,nfxx,hx,gx,gxx,eta)

save;