% Solution to Ch. 3 exercise 1 (first part)

% define vectors of different values of beta amnd gamma

npt=41 ;                         % number of points
beta=(0.9:0.2/(npt-1):1.1)' ;    % beta grid
gamm=(0.0:20.0/(npt-1):20.0)' ;  % gamma grid
gamm(1)=0.01 ;                   % first gamma > 0

% set the parameters of the law of motion of consumption growth

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

% choose sample size of simulation

ss=101 ;

% initialize matrices to store equity premia and risk free rates

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

% load the quadrature data and choose the relevant submatrix
% of data for the number of quadrature points

nqp=8 ; % number of quadrature points to use ( must be <= 10 )

load ghquad.dat ;
ymat=ghquad(0.5*(nqp-1)*nqp+1:0.5*nqp*(nqp+1),2) ;
wmat=ghquad(0.5*(nqp-1)*nqp+1:0.5*nqp*(nqp+1),3) ;

% ymat is the vector of abscissas
% wmat is the vector of weights

% Now transform the abscissas given the law of motion

ymat=mu+ymat*sig ;

% generate the transition matrix

pim=gettrans(ymat,wmat,mu,rho,sig) ; % calls external m file

% simulate discrete state space process

ix=dspace(pim,ss) ; % calls external m file to get indexes of
                    % the simulated state space indexes
x=ymat(ix,1) ;  % simulated values of consumption growth

% loop over the values of beta and gamma

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

  % get solution for price-dividend ratio

  pdr=getprice(ymat,pim,b,g) ; % calls external m file
  v=pdr(ix,1) ;  % simulated time series of price dividend ratios

  % get solution for price of risk free asset

  qp=getrfree(ymat,pim,b,g) ; % calls external m file
  q=qp(ix,1) ;   % simulated time series of risk free prices

  % 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 ;

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

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