function[r,dr_dz,d2r_dz2] = derret(zSS,pstar); % % Compute derivatives for return of the % Burnside, Eichenbaum and Rebelo model % % nz = length(zSS); np = length(pstar); A=zSS(2); g=zSS(3); k=zSS(4); n1=zSS(5); n1p1=zSS(6); w=zSS(7); kp1=zSS(8); global TE alpha=pstar(1,1); beta=pstar(2,1); delta=pstar(3,1); muA=pstar(4,1); muG=pstar(5,1); rhoA=pstar(6,1); rhoG=pstar(7,1); h1=pstar(8,1); v=pstar(9,1); xi=pstar(10,1); gamma=pstar(11,1); % a0 = -v*log(TE); a1 = -v*log(TE-xi-h1*exp(w)); z = zSS; dr_dzp = -9999*ones(nz,1); d2r_dzp2 = -9999*ones(nz,nz); % % % exc = exp(A)*exp(k)^(1-alpha)*h1*(exp(w)*exp(n1))^alpha - ... exp(g) - gamma*exp(kp1) + (1-delta)*exp(k); r = log(exc) - a1*exp(n1) - a0*(1-exp(n1)); % % derivatives computed by Maple % dr_dzp(1,1)=0; dr_dzp(2,1)=exp(A)*... exp(k)^(1-alpha)*h1*(exp(n1)*exp(w))^alpha/(exp(A)*... exp(k)^(1-alpha)*h1*(exp(n1)*exp(w))^alpha-gamma*... exp(kp1)+(1-delta)*exp(k)-exp(g)); dr_dzp(3,1)=-exp(g)/(exp(A)*exp(k)^(1-alpha)*... h1*(exp(n1)*exp(w))^alpha-gamma*exp(kp1)+(1-delta)*... exp(k)-exp(g)); dr_dzp(4,1)=(exp(A)*... exp(k)^(1-alpha)*(1-alpha)*h1*(exp(n1)*... exp(w))^alpha+(1-delta)*exp(k))/(exp(A)*... exp(k)^(1-alpha)*h1*(exp(n1)*exp(w))^alpha-gamma*... exp(kp1)+(1-delta)*exp(k)-exp(g)); dr_dzp(5,1)=exp(A)*... exp(k)^(1-alpha)*h1*(exp(n1)*exp(w))^alpha*... alpha/(exp(A)*exp(k)^(1-alpha)*h1*... (exp(n1)*exp(w))^alpha-gamma*exp(kp1)+(1-delta)*... exp(k)-exp(g))+v*log(TE-xi-exp(w)*... h1)*exp(n1)-v*log(TE)*exp(n1); dr_dzp(6,1)=0; dr_dzp(7,1)=exp(A)*... exp(k)^(1-alpha)*h1*(exp(n1)*exp(w))^alpha*... alpha/(exp(A)*exp(k)^(1-alpha)*h1*... (exp(n1)*exp(w))^alpha-gamma*exp(kp1)+(1-delta)*... exp(k)-exp(g))-v*exp(w)*h1/(TE-xi-exp(w)*... h1)*exp(n1); dr_dzp(8,1)=-gamma*exp(kp1)/(exp(A)*... exp(k)^(1-alpha)*h1*(exp(n1)*exp(w))^alpha-gamma*... exp(kp1)+(1-delta)*exp(k)-exp(g)); d2r_dzp2(6, 8)=0; d2r_dzp2(3, 4)=(exp(A)*... exp(k)^(1-alpha)*(1-alpha)*h1*(exp(n1)*... exp(w))^alpha+(1-delta)*exp(k))/(exp(A)*... exp(k)^(1-alpha)*h1*(exp(n1)*exp(w))^alpha-gamma*... exp(kp1)+(1-delta)*exp(k)-exp(g))^2*... exp(g); d2r_dzp2(4, 3)=(exp(A)*exp(k)^(1-alpha)*... (1-alpha)*h1*(exp(n1)*exp(w))^alpha+(1-delta)*... exp(k))/(exp(A)*exp(k)^(1-alpha)*... h1*(exp(n1)*exp(w))^alpha-gamma*exp(kp1)+(1-delta)*... exp(k)-exp(g))^2*exp(g); d2r_dzp2(5, 6)=0; d2r_dzp2(2, 7)=exp(A)*... exp(k)^(1-alpha)*h1*(exp(n1)*exp(w))^alpha*... alpha/(exp(A)*exp(k)^(1-alpha)*h1*... (exp(n1)*exp(w))^alpha-gamma*exp(kp1)+(1-delta)*... exp(k)-exp(g))-exp(A)^2*(exp(k)^(1-alpha))^2*... h1^2*((exp(n1)*exp(w))^alpha)^2*alpha/(exp(A)*... exp(k)^(1-alpha)*h1*(exp(n1)*exp(w))^alpha-gamma*... exp(kp1)+(1-delta)*exp(k)-exp(g))^2; d2r_dzp2(7, 1)=0; d2r_dzp2(8, 7)=gamma*... exp(kp1)/(exp(A)*exp(k)^(1-alpha)*... h1*(exp(n1)*exp(w))^alpha-gamma*exp(kp1)+(1-delta)*... exp(k)-exp(g))^2*exp(A)*exp(k)^(1-alpha)*... h1*(exp(n1)*exp(w))^alpha*alpha; d2r_dzp2(8, 8)=-gamma*exp(kp1)/(exp(A)*... exp(k)^(1-alpha)*h1*(exp(n1)*exp(w))^alpha-gamma*... exp(kp1)+(1-delta)*exp(k)-exp(g))-gamma^2*... exp(kp1)^2/(exp(A)*exp(k)^(1-alpha)*... h1*(exp(n1)*exp(w))^alpha-gamma*exp(kp1)+(1-delta)*... exp(k)-exp(g))^2; d2r_dzp2(8, 4)=gamma*... exp(kp1)/(exp(A)*exp(k)^(1-alpha)*... h1*(exp(n1)*exp(w))^alpha-gamma*exp(kp1)+(1-delta)*... exp(k)-exp(g))^2*(exp(A)*exp(k)^(1-alpha)*... (1-alpha)*h1*(exp(n1)*exp(w))^alpha+(1-delta)*... exp(k)); d2r_dzp2(1, 6)=0; d2r_dzp2(4, 6)=0; d2r_dzp2(5, 7)=exp(A)*... exp(k)^(1-alpha)*h1*(exp(n1)*exp(w))^alpha*... alpha^2/(exp(A)*exp(k)^(1-alpha)*... h1*(exp(n1)*exp(w))^alpha-gamma*exp(kp1)+(1-delta)*... exp(k)-exp(g))-exp(A)^2*(exp(k)^(1-alpha))^2*... h1^2*((exp(n1)*exp(w))^alpha)^2*alpha^2/(exp(A)*... exp(k)^(1-alpha)*h1*(exp(n1)*exp(w))^alpha-gamma*... exp(kp1)+(1-delta)*exp(k)-exp(g))^2-v*... exp(w)*h1/(TE-xi-exp(w)*h1)*exp(n1); d2r_dzp2(7, 2)=exp(A)*... exp(k)^(1-alpha)*h1*(exp(n1)*exp(w))^alpha*... alpha/(exp(A)*exp(k)^(1-alpha)*h1*... (exp(n1)*exp(w))^alpha-gamma*exp(kp1)+(1-delta)*... exp(k)-exp(g))-exp(A)^2*(exp(k)^(1-alpha))^2*... h1^2*((exp(n1)*exp(w))^alpha)^2*alpha/(exp(A)*... exp(k)^(1-alpha)*h1*(exp(n1)*exp(w))^alpha-gamma*... exp(kp1)+(1-delta)*exp(k)-exp(g))^2; d2r_dzp2(3, 1)=0; d2r_dzp2(1, 7)=0; d2r_dzp2(4, 7)=exp(A)*... exp(k)^(1-alpha)*(1-alpha)*h1*(exp(n1)*... exp(w))^alpha*alpha/(exp(A)*exp(k)^(1-alpha)*... h1*(exp(n1)*exp(w))^alpha-gamma*exp(kp1)+(1-delta)*... exp(k)-exp(g))-exp(A)*exp(k)^(1-alpha)*... h1*(exp(n1)*exp(w))^alpha*alpha/(exp(A)*... exp(k)^(1-alpha)*h1*(exp(n1)*exp(w))^alpha-gamma*... exp(kp1)+(1-delta)*exp(k)-exp(g))^2*... (exp(A)*exp(k)^(1-alpha)*(1-alpha)*... h1*(exp(n1)*exp(w))^alpha+(1-delta)*... exp(k)); d2r_dzp2(2, 3)=exp(g)/(exp(A)*... exp(k)^(1-alpha)*h1*(exp(n1)*exp(w))^alpha-gamma*... exp(kp1)+(1-delta)*exp(k)-exp(g))^2*... exp(A)*exp(k)^(1-alpha)*h1*(exp(n1)*... exp(w))^alpha; d2r_dzp2(5, 4)=exp(A)*... exp(k)^(1-alpha)*(1-alpha)*h1*(exp(n1)*... exp(w))^alpha*alpha/(exp(A)*exp(k)^(1-alpha)*... h1*(exp(n1)*exp(w))^alpha-gamma*exp(kp1)+(1-delta)*... exp(k)-exp(g))-exp(A)*exp(k)^(1-alpha)*... h1*(exp(n1)*exp(w))^alpha*alpha/(exp(A)*... exp(k)^(1-alpha)*h1*(exp(n1)*exp(w))^alpha-gamma*... exp(kp1)+(1-delta)*exp(k)-exp(g))^2*... (exp(A)*exp(k)^(1-alpha)*(1-alpha)*... h1*(exp(n1)*exp(w))^alpha+(1-delta)*... exp(k)); d2r_dzp2(1, 2)=0; d2r_dzp2(5, 1)=0; d2r_dzp2(7, 8)=gamma*exp(kp1)/(exp(A)*... exp(k)^(1-alpha)*h1*(exp(n1)*exp(w))^alpha-gamma*... exp(kp1)+(1-delta)*exp(k)-exp(g))^2*... exp(A)*exp(k)^(1-alpha)*h1*(exp(n1)*... exp(w))^alpha*alpha; d2r_dzp2(3, 2)=exp(g)/(exp(A)*... exp(k)^(1-alpha)*h1*(exp(n1)*exp(w))^alpha-gamma*... exp(kp1)+(1-delta)*exp(k)-exp(g))^2*... exp(A)*exp(k)^(1-alpha)*h1*(exp(n1)*... exp(w))^alpha; d2r_dzp2(3, 5)=exp(A)*... exp(k)^(1-alpha)*h1*(exp(n1)*exp(w))^alpha*... alpha/(exp(A)*exp(k)^(1-alpha)*h1*... (exp(n1)*exp(w))^alpha-gamma*exp(kp1)+(1-delta)*... exp(k)-exp(g))^2*exp(g); d2r_dzp2(1, 4)=0; d2r_dzp2(6, 5)=0; d2r_dzp2(2, 8)=gamma*... exp(kp1)/(exp(A)*exp(k)^(1-alpha)*... h1*(exp(n1)*exp(w))^alpha-gamma*exp(kp1)+(1-delta)*... exp(k)-exp(g))^2*exp(A)*exp(k)^(1-alpha)*... h1*(exp(n1)*exp(w))^alpha; d2r_dzp2(4, 4)=(exp(A)*exp(k)^(1-alpha)*... (1-alpha)^2*h1*(exp(n1)*exp(w))^alpha+(1-delta)*... exp(k))/(exp(A)*exp(k)^(1-alpha)*... h1*(exp(n1)*exp(w))^alpha-gamma*exp(kp1)+(1-delta)*... exp(k)-exp(g))-(exp(A)*exp(k)^(1-alpha)*... (1-alpha)*h1*(exp(n1)*exp(w))^alpha+(1-delta)*... exp(k))^2/(exp(A)*exp(k)^(1-alpha)*... h1*(exp(n1)*exp(w))^alpha-gamma*exp(kp1)+(1-delta)*... exp(k)-exp(g))^2; d2r_dzp2(7, 5)=exp(A)*exp(k)^(1-alpha)*... h1*(exp(n1)*exp(w))^alpha*alpha^2/(exp(A)*... exp(k)^(1-alpha)*h1*(exp(n1)*exp(w))^alpha-gamma*... exp(kp1)+(1-delta)*exp(k)-exp(g))-exp(A)^2*... (exp(k)^(1-alpha))^2*h1^2*((exp(n1)*... exp(w))^alpha)^2*alpha^2/(exp(A)*... exp(k)^(1-alpha)*h1*(exp(n1)*exp(w))^alpha-gamma*... exp(kp1)+(1-delta)*exp(k)-exp(g))^2-v*... exp(w)*h1/(TE-xi-exp(w)*h1)*exp(n1); d2r_dzp2(5, 2)=exp(A)*... exp(k)^(1-alpha)*h1*(exp(n1)*exp(w))^alpha*... alpha/(exp(A)*exp(k)^(1-alpha)*h1*... (exp(n1)*exp(w))^alpha-gamma*exp(kp1)+(1-delta)*... exp(k)-exp(g))-exp(A)^2*(exp(k)^(1-alpha))^2*... h1^2*((exp(n1)*exp(w))^alpha)^2*alpha/(exp(A)*... exp(k)^(1-alpha)*h1*(exp(n1)*exp(w))^alpha-gamma*... exp(kp1)+(1-delta)*exp(k)-exp(g))^2; d2r_dzp2(8, 1)=0; d2r_dzp2(8, 5)=gamma*... exp(kp1)/(exp(A)*exp(k)^(1-alpha)*... h1*(exp(n1)*exp(w))^alpha-gamma*exp(kp1)+(1-delta)*... exp(k)-exp(g))^2*exp(A)*exp(k)^(1-alpha)*... h1*(exp(n1)*exp(w))^alpha*alpha; d2r_dzp2(6, 6)=0; d2r_dzp2(4, 8)=gamma*... exp(kp1)/(exp(A)*exp(k)^(1-alpha)*... h1*(exp(n1)*exp(w))^alpha-gamma*exp(kp1)+(1-delta)*... exp(k)-exp(g))^2*(exp(A)*exp(k)^(1-alpha)*... (1-alpha)*h1*(exp(n1)*exp(w))^alpha+(1-delta)*... exp(k)); d2r_dzp2(5, 8)=gamma*... exp(kp1)/(exp(A)*exp(k)^(1-alpha)*... h1*(exp(n1)*exp(w))^alpha-gamma*exp(kp1)+(1-delta)*... exp(k)-exp(g))^2*exp(A)*exp(k)^(1-alpha)*... h1*(exp(n1)*exp(w))^alpha*alpha; d2r_dzp2(8, 2)=gamma*... exp(kp1)/(exp(A)*exp(k)^(1-alpha)*... h1*(exp(n1)*exp(w))^alpha-gamma*exp(kp1)+(1-delta)*... exp(k)-exp(g))^2*exp(A)*exp(k)^(1-alpha)*... h1*(exp(n1)*exp(w))^alpha; d2r_dzp2(4, 1)=0; d2r_dzp2(7, 4)=exp(A)*... exp(k)^(1-alpha)*(1-alpha)*h1*(exp(n1)*... exp(w))^alpha*alpha/(exp(A)*exp(k)^(1-alpha)*... h1*(exp(n1)*exp(w))^alpha-gamma*exp(kp1)+(1-delta)*... exp(k)-exp(g))-exp(A)*exp(k)^(1-alpha)*... h1*(exp(n1)*exp(w))^alpha*alpha/(exp(A)*... exp(k)^(1-alpha)*h1*(exp(n1)*exp(w))^alpha-gamma*... exp(kp1)+(1-delta)*exp(k)-exp(g))^2*... (exp(A)*exp(k)^(1-alpha)*(1-alpha)*... h1*(exp(n1)*exp(w))^alpha+(1-delta)*... exp(k)); d2r_dzp2(6, 7)=0; d2r_dzp2(1, 1)=0; d2r_dzp2(2, 4)=exp(A)*... exp(k)^(1-alpha)*(1-alpha)*h1*(exp(n1)*... exp(w))^alpha/(exp(A)*exp(k)^(1-alpha)*... h1*(exp(n1)*exp(w))^alpha-gamma*exp(kp1)+(1-delta)*... exp(k)-exp(g))-(exp(A)*exp(k)^(1-alpha)*... (1-alpha)*h1*(exp(n1)*exp(w))^alpha+(1-delta)*... exp(k))/(exp(A)*exp(k)^(1-alpha)*... h1*(exp(n1)*exp(w))^alpha-gamma*exp(kp1)+(1-delta)*... exp(k)-exp(g))^2*exp(A)*exp(k)^(1-alpha)*... h1*(exp(n1)*exp(w))^alpha; d2r_dzp2(6, 1)=0; d2r_dzp2(3, 3)=-exp(g)/(exp(A)*... exp(k)^(1-alpha)*h1*(exp(n1)*exp(w))^alpha-gamma*... exp(kp1)+(1-delta)*exp(k)-exp(g))-exp(g)^2/(exp(A)*... exp(k)^(1-alpha)*h1*(exp(n1)*exp(w))^alpha-gamma*... exp(kp1)+(1-delta)*exp(k)-exp(g))^2; d2r_dzp2(4, 2)=exp(A)*... exp(k)^(1-alpha)*(1-alpha)*h1*(exp(n1)*... exp(w))^alpha/(exp(A)*exp(k)^(1-alpha)*... h1*(exp(n1)*exp(w))^alpha-gamma*exp(kp1)+(1-delta)*... exp(k)-exp(g))-(exp(A)*exp(k)^(1-alpha)*... (1-alpha)*h1*(exp(n1)*exp(w))^alpha+(1-delta)*... exp(k))/(exp(A)*exp(k)^(1-alpha)*... h1*(exp(n1)*exp(w))^alpha-gamma*exp(kp1)+(1-delta)*... exp(k)-exp(g))^2*exp(A)*exp(k)^(1-alpha)*... h1*(exp(n1)*exp(w))^alpha; d2r_dzp2(5, 5)=exp(A)*... exp(k)^(1-alpha)*h1*(exp(n1)*exp(w))^alpha*... alpha^2/(exp(A)*exp(k)^(1-alpha)*... h1*(exp(n1)*exp(w))^alpha-gamma*exp(kp1)+(1-delta)*... exp(k)-exp(g))-exp(A)^2*(exp(k)^(1-alpha))^2*... h1^2*((exp(n1)*exp(w))^alpha)^2*alpha^2/(exp(A)*... exp(k)^(1-alpha)*h1*(exp(n1)*exp(w))^alpha-gamma*... exp(kp1)+(1-delta)*exp(k)-exp(g))^2+v*... log(TE-xi-exp(w)*h1)*exp(n1)-v*log(TE)*... exp(n1); d2r_dzp2(3, 6)=0; d2r_dzp2(6, 2)=0; d2r_dzp2(7, 3)=exp(A)*... exp(k)^(1-alpha)*h1*(exp(n1)*exp(w))^alpha*... alpha/(exp(A)*exp(k)^(1-alpha)*h1*... (exp(n1)*exp(w))^alpha-gamma*exp(kp1)+(1-delta)*... exp(k)-exp(g))^2*exp(g); d2r_dzp2(8, 3)=-gamma*... exp(kp1)/(exp(A)*exp(k)^(1-alpha)*... h1*(exp(n1)*exp(w))^alpha-gamma*exp(kp1)+(1-delta)*... exp(k)-exp(g))^2*exp(g); d2r_dzp2(8, 6)=0; d2r_dzp2(3, 7)=exp(A)*... exp(k)^(1-alpha)*h1*(exp(n1)*exp(w))^alpha*... alpha/(exp(A)*exp(k)^(1-alpha)*h1*... (exp(n1)*exp(w))^alpha-gamma*exp(kp1)+(1-delta)*... exp(k)-exp(g))^2*exp(g); d2r_dzp2(1, 5)=0; d2r_dzp2(1, 8)=0; d2r_dzp2(1, 3)=0; d2r_dzp2(7, 6)=0; d2r_dzp2(2, 5)=exp(A)*... exp(k)^(1-alpha)*h1*(exp(n1)*exp(w))^alpha*... alpha/(exp(A)*exp(k)^(1-alpha)*h1*... (exp(n1)*exp(w))^alpha-gamma*exp(kp1)+(1-delta)*... exp(k)-exp(g))-exp(A)^2*(exp(k)^(1-alpha))^2*... h1^2*((exp(n1)*exp(w))^alpha)^2*alpha/(exp(A)*... exp(k)^(1-alpha)*h1*(exp(n1)*exp(w))^alpha-gamma*... exp(kp1)+(1-delta)*exp(k)-exp(g))^2; d2r_dzp2(4, 5)=exp(A)*... exp(k)^(1-alpha)*(1-alpha)*h1*(exp(n1)*... exp(w))^alpha*alpha/(exp(A)*exp(k)^(1-alpha)*... h1*(exp(n1)*exp(w))^alpha-gamma*exp(kp1)+(1-delta)*... exp(k)-exp(g))-exp(A)*exp(k)^(1-alpha)*... h1*(exp(n1)*exp(w))^alpha*alpha/(exp(A)*... exp(k)^(1-alpha)*h1*(exp(n1)*exp(w))^alpha-gamma*... exp(kp1)+(1-delta)*exp(k)-exp(g))^2*... (exp(A)*exp(k)^(1-alpha)*(1-alpha)*... h1*(exp(n1)*exp(w))^alpha+(1-delta)*... exp(k)); d2r_dzp2(6, 3)=0; d2r_dzp2(2, 1)=0; d2r_dzp2(7, 7)=exp(A)*... exp(k)^(1-alpha)*h1*(exp(n1)*exp(w))^alpha*... alpha^2/(exp(A)*exp(k)^(1-alpha)*... h1*(exp(n1)*exp(w))^alpha-gamma*exp(kp1)+(1-delta)*... exp(k)-exp(g))-exp(A)^2*(exp(k)^(1-alpha))^2*... h1^2*((exp(n1)*exp(w))^alpha)^2*alpha^2/(exp(A)*... exp(k)^(1-alpha)*h1*(exp(n1)*exp(w))^alpha-gamma*... exp(kp1)+(1-delta)*exp(k)-exp(g))^2-v*... exp(w)*h1/(TE-xi-exp(w)*h1)*exp(n1)-v*... exp(w)^2*h1^2/(TE-xi-exp(w)*h1)^2*... exp(n1); d2r_dzp2(5, 3)=exp(A)*... exp(k)^(1-alpha)*h1*(exp(n1)*exp(w))^alpha*... alpha/(exp(A)*exp(k)^(1-alpha)*h1*... (exp(n1)*exp(w))^alpha-gamma*exp(kp1)+(1-delta)*... exp(k)-exp(g))^2*exp(g); d2r_dzp2(2, 2)=exp(A)*... exp(k)^(1-alpha)*h1*(exp(n1)*exp(w))^alpha/(exp(A)*... exp(k)^(1-alpha)*h1*(exp(n1)*exp(w))^alpha-gamma*... exp(kp1)+(1-delta)*exp(k)-exp(g))-exp(A)^2*... (exp(k)^(1-alpha))^2*h1^2*((exp(n1)*... exp(w))^alpha)^2/(exp(A)*exp(k)^(1-alpha)*... h1*(exp(n1)*exp(w))^alpha-gamma*exp(kp1)+(1-delta)*... exp(k)-exp(g))^2; d2r_dzp2(6, 4)=0; d2r_dzp2(3, 8)=-gamma*exp(kp1)/(exp(A)*... exp(k)^(1-alpha)*h1*(exp(n1)*exp(w))^alpha-gamma*... exp(kp1)+(1-delta)*exp(k)-exp(g))^2*... exp(g); d2r_dzp2(2, 6)=0; % % % dr_dz = dr_dzp; d2r_dz2 = d2r_dzp2;