function v=getprice(ymat,pim,bet,gam) ;
% function to solve the euler equation
                   % ymat is the vector of quadrature points
                   % pim is the transition matrix
                   % bet is beta
                   % gam is gamma

[n,n0]=size(pim) ; % determines the number of points in state space = n

mrs=(bet*exp(ymat(:,ones(n,1))*(1-gam)))' ;
                   % computes an nxn matrix whose ji element is
                   % beta*exp[(1-gamma)*y(i)]. So it is the value of the
                   % variable beta*exp[(1-gamma)*x(t+1)] if state j
                   % occured at time t, and i occured at time t+1 -
                   % notice it's independent of j.

pm=pim.*mrs ;      % creates the element product of pim and mrs

% the Euler equation is  v = pm*(1+v) because pm*(1+v) creates
% sum(i) pm(j,i)*[1+v(i)].  So the Euler equation can be rewritten as
% (I-pm)*v = pm*1 for a unit vector 1.  The next statement solves this
% equation for v

v=(eye(n)-pm)\(pm*ones(n,1)) ;
