function[r,dr_dz,d2r_dz2] = derret(zSS,pstar)
%
%   Compute derivatives for return of the 
%           Hansen (1985) model
%
%
nz = length(zSS);
lambda=zSS(2); kss=zSS(3); iss=zSS(4); hss=zSS(5); 

theta = pstar(1,1); 
delta = pstar(2,1);  
beta  = pstar(3,1);   
A     = pstar(4,1);
rho   = pstar(5,1);   
%
%
%
dr_dzp = -9999*ones(nz,1);
d2r_dz2 = -9999*ones(nz,nz);
%
%
%
css = lambda*(kss^theta)*hss^(1-theta) - iss;
r   = log(css) + A*log(1-hss);
%
%    first derivatives
%
dr_dz(1,1) = 0;
dr_dz(2,1) = (kss^theta*hss^(1-theta))/css;
dr_dz(3,1) = (theta*lambda*kss^(theta-1)*hss^(1-theta))/css;
dr_dz(4,1) = -1/css;
dr_dz(5,1) = ((1-theta)*lambda*kss^(theta)*hss^(-theta))/css - A/(1-hss);
%
%   compute second derivatives numerically
%
eps = 0.001;
for i = 1:nz
    zh = zSS;
    zl = zSS;
    zh(i) = zh(i) + eps;
    zl(i) = zl(i) - eps; 
    [ r_high, dr_dz_high ] = rfunc(zh,pstar);
    [ r_low, dr_dz_low ]  = rfunc(zl,pstar);
    
    dr_dzp(i) = (r_high - r_low)/(2*eps);
    d2r_dz2(:,i) = (dr_dz_high - dr_dz_low)/(2*eps);
end