% Program that solves Ch. 2 Exercise 2

% define the law of motion parameters

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

% set model parameters

b=0.99 ;
g=2.0 ;

% set range of the x variable

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

% get the linearized solutions

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

ql=exp(qt-g*rho*(x-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)+0.5*g*g*sigma*sigma) ;
v=zeros(101,1) ;

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

while ((i <= 10000) & (meps > 0.000001)) ;
 eps=betai*exp(a*(mu*(i-rho*(1-rhoi)/(1-rho))+ ...
   (rho*(1-rhoi)/(1-rho))*x)+0.5*a*a*sigma*sigma*(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 ;
i

plot(x,[v vl]) ;
xlabel('x')
ylabel('V')
title('True (yellow) Log-lin (magenta)')
pause

plot(x,v./vl) ;
xlabel('x')
title('Ratio of True to Log-lin')

