% function [cp, w, p,a,b] = getwp(e,sd,r,nyrs,n);
e = .05;
sd = .09;
r = .02;
nyrs = 10;
n = 100;


%   function [cp, w, p,a,b] = getwp(e,sd,r,nyrs,n);
%   gets cumulative probabilities, wealth levels and state prices 
%  inputs:
%        e:  annual expected return (not in logs) (e.g. 0.10)
%       sd:  annual standard deviation of return (not in logs) (e.g. 0.18)
%        r:  annual risk-free rate (e.g. 0.04)
%     nyrs: number of years for forecast
%        n:  number of discrete outcomes
%  outputs:
%       cp:  cumulative probabilities
%        w:  ending wealth levels
%        p:  state prices
%        a:  coefficient in log(p) = a + b*log(w)
%        b:  coefficient in log(p) = a + b*log(w)

% compute riskless value relative
   rf =  1+r;

% compute moments of lognormal distribution (a is mean, b is sd)
     u = 1 + e;
     v = sd^2;
     b = sqrt(log((v/(u^2)) + 1)); 
     a = 0.5 *log((u^2)/exp(b^2));

% assume that wealth relative is lognormal with mean lwbar and standard deviation lwsd
    lwbar = a*nyrs;
    lwsd  = b*sqrt(nyrs);

% get discrete points from normal curve
   dev =  discnorm(n);

% find wealth values
   dd = lwbar + dev*lwsd;
   w = exp(dd);

% make vector of cumulative probabilities of w
   cp = 1/(2*n):1/n:1;

% find state prices such that log(p) = a + b*log(w)
%   and sum(pw) = 1
  % save w, rf and nyrs
       save temp w rf nyrs;
  % get value of b
       b = fzero('getb',0);
      % note: getb is as follows:
      %     function zval = getb(b)
      %     load temp;  % gets w, rf and nyrs
      %     zval = sum(w.^(b+1)) - (rf^nyrs)*sum(w.^b);
  % get value of a
      a = log(1/sum(w.^(b+1)));
  % compute prices
      p = exp(a)*(w.^b);


