% Program that solves Ch. 2 Exercise 1

% define vectors of different values of beta amnd gamma

npt=41 ;
beta=(0.9:0.2/(npt-1):1.1)' ;
gamm=(0.0:20.0/(npt-1):20.0)' ;
gamm(1)=0.01 ;

% define the law of motion parameters

mu=0.0179 ;
rho=-0.139 ;
sigma=0.0348 ;

% choose sample size of simulation-this is short for illustration

ss=101 ;

% initialize matrices to store equity premia and risk free rates

bigepr=zeros(npt,npt) ;
bigrfr=zeros(npt,npt) ;

% simulate consumption (dividend) growth

x=normar1(mu,rho,sigma,ss) ;

% loop over the values of beta and gamma

for i=1:npt ;
 for j=1:npt ;
  b=beta(i,1) ;
  g=gamm(j,1) ;


% solve steady state conditions
% a trap must be built in to catch inadmissable pairs of beta and gamma

  if (b*exp((1-g)*mu)) < 1,
   vt=log(b)+(1-g)*mu-log(1-b*exp((1-g)*mu)) ;
   qt=log(b)-g*mu ;
   m=exp(log(b)+(1-g)*mu) ;

   % compute simulated prices

   v=exp(vt+(1-g)*rho*(x-mu)/(1-m*rho)) ;
   q=exp(qt-g*rho*(x-mu)) ;

   % compute rates of return

   re=(v(2:ss,1)+1).*exp(x(2:ss,1))./v(1:ss-1,1)-1 ;
   rf=ones(ss-1,1)./q(1:ss-1,1)-1 ;

   % compute and store the averages

   rfr=mean(rf) ;
   epr=mean(re)-rfr ;
   bigrfr(i,j)=100*rfr ;
   bigepr(i,j)=100*epr ;

  end ;

 end ;
end ;

surf(beta,gamm,bigepr') ;
title('The Equity Premium') ;
zlabel('%')
ylabel('gamma')
xlabel('beta') ;

surf(beta,gamm,bigrfr') ;
title('The Risk Free Rate') ;
zlabel('%')
ylabel('gamma')
xlabel('beta') ;
