function [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,approx); %[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,approx); % Computes analytical first and second (if approx=2) derivatives of the function f(yp,y,xp,x) with respect to x, y, xp, and yp. For documentation, see the paper ``Solving Dynamic General Equilibrium Models Using a Second-Order Approximation to the Policy Function,'' by Stephanie Schmitt-Grohe and Martin Uribe, 2001). % %Inputs: f, x, y, xp, yp, approx % %Output: Analytical first and second derivatives of f. % %If approx is set at a value different from 2, the program delivers the first derivatives of f and sets second derivatives at zero. If approx equals 2, the program returns first and second derivatives of f. The default value of approx is 2. %Note: This program requires MATLAB's Symbolic Math Toolbox % %(c) Stephanie Schmitt-Grohe and Martin Uribe %Date July 17, 2001 if nargin==5 approx=2; end nx = size(x,2); ny = size(y,2); nxp = size(xp,2); nyp = size(yp,2); n = size(f,1); %Compute the first and second derivatives of f fx = jacobian(f,x); fxp = jacobian(f,xp); fy = jacobian(f,y); fyp = jacobian(f,yp); if approx==2 fypyp = reshape(jacobian(fyp(:),yp),n,nyp,nyp); fypy = reshape(jacobian(fyp(:),y),n,nyp,ny); fypxp = reshape(jacobian(fyp(:),xp),n,nyp,nxp); fypx = reshape(jacobian(fyp(:),x),n,nyp,nx); fyyp = reshape(jacobian(fy(:),yp),n,ny,nyp); fyy = reshape(jacobian(fy(:),y),n,ny,ny); fyxp = reshape(jacobian(fy(:),xp),n,ny,nxp); fyx = reshape(jacobian(fy(:),x),n,ny,nx); fxpyp = reshape(jacobian(fxp(:),yp),n,nxp,nyp); fxpy = reshape(jacobian(fxp(:),y),n,nxp,ny); fxpxp = reshape(jacobian(fxp(:),xp),n,nxp,nxp); fxpx = reshape(jacobian(fxp(:),x),n,nxp,nx); fxyp = reshape(jacobian(fx(:),yp),n,nx,nyp); fxy = reshape(jacobian(fx(:),y),n,nx,ny); fxxp = reshape(jacobian(fx(:),xp),n,nx,nxp); fxx = reshape(jacobian(fx(:),x),n,nx,nx); else fypyp=0; fypy=0; fypxp=0; fypx=0; fyyp=0; fyy=0; fyxp=0; fyx=0; fxpyp=0; fxpy=0; fxpxp=0; fxpx=0; fxyp=0; fxy=0; fxxp=0; fxx=0; end