% Solution to Ch 3 exercise 3

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

mu=0.0179 ;
rho=-0.139 ;
sig=0.0348 ;
sigx=sqrt(sig*sig/(1-rho*rho)) ;

% choose sample size of simulation

ss=1001 ;

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

load ghquad.dat ;

% put in some output labels

'    No.Pts.          Mean             Std. Dev.'
'               R.F.      E.P.       R.F.       E.P.'

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

 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*sigx ;

 % generate the transition matrix

 pim=gettran2(ymat,wmat,mu,rho,sig,sigx) ;
 pu=ones(nqp,1)./nqp ;
 md=1 ;
 while md > 1e-12 ;
  pup=pim'*pu ;
  md=max(abs(pup-pu)) ;
  pu=pup ;
 end ;

 % choose values of beta and gamma

 b=0.99 ;
 g=2.0 ;

 % get solution for price dividend ratio

 pdr=getprice(ymat,pim,b,g) ;
 eqr=(1+pdr(:,ones(nqp,1))').*exp(ymat(:,ones(nqp,1))')./ ...
     pdr(:,ones(nqp,1))-1 ;

 % get solution for price of risk free asset

 qp=getrfree(ymat,pim,b,g) ;
 rfr=ones(nqp,1)./qp-1 ;
 rfr=rfr(:,ones(nqp,1)) ;

 % compute and print the moments

 meanrfr=sum((pim.*rfr)')*pu ;
 meanepr=sum((pim.*(eqr-rfr))')*pu ;
 stdrfr=sqrt(sum((pim.*(rfr-meanrfr).*(rfr-meanrfr))')*pu) ;
 stdepr=sqrt(sum((pim.*(eqr-rfr-meanepr).*(eqr-rfr-meanepr))')*pu) ;

 [nqp meanrfr meanepr stdrfr stdepr]

end ;
