function [y,theta,k, cumprob,xc] = dgamma(x,m,s,pc);
%y = dgamma(x,m,s) evaluates at x the 
%probability density function of the gamma 
%distribution with mean m and standard deviation  s.
%[y,theta,k] = dgamma(x,m,s) produces the 
%parameters defining the probability density 
%x^(k-1)*exp(-x/theta)/gamma(k)/(theta^k);
%%[y,theta,k, cumprob,xc] = dgamma(x,m,s) produces the value cumprob  representing the probability that a random variable with gamma  distribution defined by the parameters theta and k is less than or equal to x. 
%[y,theta,k, cumprob,xc] = dgamma(x,m,s,pc) takes as input the value pc and produces the output xc such that the probability that a random variable with distribution gamma and parameters theta k be less than xc is pc. So, for example, if pc=0.5, then xc is themedian of the gamma  distribution
%if x and m are negative the output returns [y,theta,k, cumprob,xc], where   [y,theta,k, cumprob,-xc] = dgamma(-x,-m,s,pc)
%Replication materials for: ``What Do Long Data Tell Us About the Inflation Hike Post COVID-19 Pandemic?''
%© Stephanie Schmitt-Grohé and Martín Uribe, December 27, 2023. 

sx = sign(x);
x = abs(x);
sm = sign(m);
m = abs(m);

theta = s.^2 ./ m;
k = m ./ theta;

y = x.^(k-1).*exp(-x./theta)./gamma(k)./(theta.^k);
y(x<0) = 0;

if nargout>3
cumprob = gamcdf(x,k,theta);
cumprob(x<0) = 0;
end

if nargin>3 
xcH = 1e16;
xcL = 0;
pc0 = inf;
while abs(pc0 - pc)>1e-8
xc = (xcH+xcL)/2;
pc0 =   gamcdf(xc,k,theta);
if pc0 > pc
xcH = xc;
else xcL = xc;
end %if ans > pc
end  %while abs(pc0 - pc)>1e-8
xc = xc.*(sm>0) + (1-xc).*(sm<0);
end; %if nargin >3

y = y .* (sx==sm);

