function [y,a,b,cumprob,xc] = dbeta(x,m,s,pc);
%y = dbeta(x,m,s) evaluates at x the 
%probability density function of the beta 
%distribution with mean m and standard deviation  s.
%[y,a,b] = dbeta(x,m,s) produces the parametes a and b defining the density functon y = x^(a-1) * (1-x)^(b-1) / beta(a,b)
%%[y,a,b,cumprob] = dbeta(x,m,s) produces the value cumprob  representing the probability that a random viarialbe with beta distribution defined by the parameters a and b is less than or equal to x. 
%[y,a,b,cumprob,xc] = dbeta(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 beta (a,b) be less than xc is pc. So, for example, if pc=0.5, then xc is themedian of the beta distribution
%Replication materials for: ``What Do Long Data Tell Us About the Permanent Component of Inflation?''
%© Stephanie Schmitt-Grohé and Martín Uribe,  2023. 


boa = (1-m) ./ m;
a = (boa ./ s.^2 ./ (1+boa).^2 -1) .* m;
b = boa .* a;

y = x.^(a-1) .* (1-x).^(b-1) ./ beta(a,b);
y(x<0|x>1) = 0;

if nargout>3
cumprob = 1./a.*x.^a.*hypergeom([a, -b+1],[1+a],x) ./ beta(a,b);
cumprob(x<0|x>1) = 0;
end

if nargin>3 

xcH = 1;
xcL = 0;
pc0 = inf;
while abs(pc0 - pc)>1e-8
xc = (xcH+xcL)/2;
pc0 =   1./a.*xc.^a.*hypergeom([a, -b+1],[1+a],xc) ./ beta(a,b);
if pc0 > pc
xcH = xc;
else xcL = xc;
end %if ans > pc
end  %while abs(pc0 - pc)>1e-8

else %if nargin>3
xc = [];
end; %if nargin >3
