% asset pricing example with serially correlated shocks



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

% calibrated to annual U.S. data from Mehra and Prescott 1985



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

zmat=ghquad(0.5*(nqp-1)*nqp+1:0.5*nqp*(nqp+1),4) ;

zmat2=[ -100 ; zmat(1:nqp-1,1) ] ;



% ymat is the vector of abscissas (points)

% wmat is the vector of weights

% zmat is the vector of upper end points of the intervals

% zmat2 is the vector of lower end points of the intervals



% Now transform the abscissas and end points given the law of motion



ymat=mu+ymat*sig ;

zmat=mu+zmat*sig ;

zmat2=mu+zmat2*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 step function extension



vstep=zeros(101,1) ;

for xi=1:101 ;

 for vi=1:nqp ;

  if (x(xi,1)<zmat(vi,1))&(x(xi,1)>=zmat2(vi,1))

   vstep(xi,1)=pdr(vi,1) ;

  end ;

 end ;

end ;



% 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 vstep vn],ymat,pdr,'cx')

title('True (blue) Step (green) Nystrom (red)')
