% Solution to Ch 3 exercise 1 (fifth part)

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

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

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

load ghquad.dat ;

nqp=8 ; % 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*sig ;

% generate the transition matrix

pim=gettrans(ymat,wmat,mu,rho,sig) ;

% choose values of beta and gamma

b=0.99 ;
g=2.0 ;

% get solution for price dividend ratio

pdr=getprice(ymat,pim,b,g) ;

% set range of the x variable

sigmax=sig/sqrt(1-rho*rho) ;
x=(-3*sigmax+mu:6*sigmax/100:3*sigmax+mu)' ;

% get the nystrom extension

vn=zeros(101,1) ;
for xi=1:101 ;
 for vi=1:nqp ;
  vn(xi,1)=vn(xi,1)+b*exp((1-g)*ymat(vi,1)).* ...
       condnorm(ymat(vi,1),x(xi,1),mu,rho,sig)./ ...
       condnorm(ymat(vi,1),mu,mu,rho,sig).*(1+pdr(vi,1)).*wmat(vi,1) ;
 end ;
end ;

% get the linearized solutions

m=exp(log(b)+(1-g)*mu) ;
vt=log(b)+(1-g)*mu-log(1-b*exp((1-g)*mu)) ;

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

% get the true solutions

q=b*exp( -g*mu-g*rho*(x-mu) ) ;
v=zeros(101,1) ;

a=1-g ;
betai=b ;
rhoi=rho ;
rho2i=rho*rho ;
i=1 ;
meps=1 ;

while ((i <= 10000) & (meps > 0.001)) ;
 eps=betai*exp(a*(mu*(i-rho*(1-rhoi)/(1-rho))+ ...
   (rho*(1-rhoi)/(1-rho))*x)+0.5*a*a*sig*sig*(i-2*rho*(1-rhoi)/ ...
   (1-rho)+rho*rho*(1-rho2i)/(1-rho*rho))/((1-rho)*(1-rho))) ;
 meps=max(abs(eps)) ;
 v=v+eps ;
 betai=b*betai ;
 rhoi=rho*rhoi ;
 rho2i=rho*rho*rho2i ;
 i=i+1 ;
end ;

plot(x,[v vl vn]) ;
xlabel('x')
ylabel('V')
title('True (yel) Loglin (mag) Nystrom (blue)')
