% example of Tauchen solution to the growth model



% define parameters



delt=0.025 ;

bet=1.03^(-0.25) ;

thet=0.4 ;

rho=0.9 ;

sigma=0.01 ;



% start off by getting the log-linear solution



% define crucial ratios in nonstochastic steady state



'K-Y Ratio'

ky=thet/(1/bet-(1-delt))

'C-Y Ratio'

cy=1-delt*ky

'C-K Ratio'

ck=cy/ky



% assume A=1 in nonstochastic steady state



kss=ky^(1/(1-thet)) ; % steady state k

css=ck*kss ;          % steady state c



mu=bet*thet/ky ;



% solve for the kappas and the alphas



k1 = (thet-1)*mu*ck-(1+1/bet) ;

k2 = 1/bet ;

a1 = (1-mu)/ky+mu*delt ;

a2 = -1/ky ;



% solve for the lambdas



'Roots of the Difference Equation'



l1=(-k1-sqrt(k1*k1-4/bet))/2 ;

l2=(-k1+sqrt(k1*k1-4/bet))/2 ;

[l1 l2]



% solve for the decision rule coefficients



'Coefficients on k(t) and a(t)'



h1=l1 ;

h2=-(a1*rho+a2)*(1/l2)/(1-rho/l2) ;

[h1 h2]



% get moments for setting up the grid for k

% do this analytically using the state space rep of the model



% | k(t+1) |   | h1    h2 | | k(t) |   |    0     |

% |        | = |          | |      | + |          |

% | a(t+1) |   | 0    rho | | a(t) |   | eps(t+1) |

%

% or

%

%  m(t+1) = M m(t) + e(t+1)

%

% which implies V(m) = V(e)+ M V(e) M' + M^2 V(e) (M^2)' + ....

% which is generated recursively below

% V(k)=V(m)(1,1)



ccrit=1 ;                              % convergence criterion

sigm=zeros(2,2) ;                      % initialize V(m)

sige=[ 0.0 0.0 ; 0.0 (sigma*sigma) ] ; % initialize V(e)

mmat=[ h1 h2 ; 0.0 rho ] ;             % initialize M

gmat=eye(2) ;                          % initialize G=M^i as I

while (ccrit > 1e-12) ;

 badd=gmat*sige*(gmat') ;              % M^i V(e) (M^i)'

 ccrit=max(max(abs(badd))') ;          % update criterion

 sigm=sigm+badd ;                      % new V(m)

 gmat=mmat*gmat ;                      % new G

end ;

sdk=sqrt(sigm(1,1)) ; % standard deviation of log(k)



% set up the grid for k



lkss=log(kss) ; % log of nssv of K for A=1 case



nkp=90 ;

kgrid=(-4*sdk:8*sdk/(nkp-1):+4*sdk)' ; % don't add the mean yet



% generate the grid for the technology shocks



nap=10 ; % number of points for the A grid



load ghquad.dat ;

agrid=ghquad(0.5*(nap-1)*nap+1:0.5*nap*(nap+1),2) ;

wmat=ghquad(0.5*(nap-1)*nap+1:0.5*nap*(nap+1),3) ;

agrid=agrid*sigma ;

pim=gettrans(agrid,wmat,0,rho,sigma) ;



% initialize the value function by obtaining the value function

% for the log-linear case



% V(k_0,A_0)=E_0 sum(t) b^t log(c_t)

%  but from the log-linear case we have the log-lin budget constraint

%  log(c_t)=log(css)+(y/c) ahat_t + [(y/c)*thet+(1-delt)*(k/c)] khat_t

%           -(k/c) khat_(t+1)

%    = log(css)+ (y/c) ahat_t + [(y/c)*thet+(1-delt)*(k/c)] khat_t

%           -(k/c) [ h1 khat_t + h2 ahat_t ]

%    = log(css)+ [(y/c)-(k/c)h2] ahat_t +

%        [(y/c)*thet+(1-delt)*(k/c)-(k/c)h1] hkat_t

%  = log(css)+ { [(y/c)*thet+(1-delt)*(k/c)-(k/c)h1] [(y/c)-(k/c)h2]  }

%     m(t) (see notation above)

%  = log(css) + phi'm(t)

% So E_0 log(c_t) = log(css)+phi M^t m(0)

% and E_0 sum(t) b^t log(c_t) = sum(t) b^t [log(css)+phi'M^t m(0)]

%  = log(css)/(1-b) + phi'inv(I-bM)*m(0)



phi1=thet/cy+(1-delt)/ck-h1/ck ;

phi2=1/cy-h2/ck ;

phi=[ phi1 phi2 ] ;

vfn=log(css)/(1-bet)+zeros(nkp,nap) ;



for ikp=1:nkp ;

 for iap=1:nap ;



  m0=[ kgrid(ikp,1) agrid(iap,1)]' ;

  vfn(ikp,iap)=vfn(ikp,iap)+phi*inv(eye(2)-bet*mmat)*m0 ;



 end ;

end ;



% initialize the decision rule



drul=zeros(nkp,nap) ;



kgrid=kgrid+lkss ; % now add the mean to the k grid



iter=0 ;

ccrit=1 ;

while ((iter <= 10000)&(ccrit > 1e-5)) ;



 vnew=zeros(nkp,nap) ; % initialize the new value function



 for ikp=1:nkp ;

  for iap=1:nap ;



   % at the current point in the state space get all possible

   % values of current consumption - make those values that are

   % negative a tiny positive number



   ccurr=exp(agrid(iap,1)+thet*kgrid(ikp,1))+(1-delt)* ...

         exp(kgrid(ikp,1))-exp(kgrid) ;

   ccurp=max([ccurr (0.000000001*ones(nkp,1)) ]')' ;



   % all possible values for this period are the logarithms of the

   % possible consumptions plus the discounted possible future

   % values



   vs=log(ccurp)+bet*( sum( (vfn.*pim(iap*ones(nkp,1),:))' )' ) ;



   % the new value function is the max of these possible values



   vnew(ikp,iap)=max(vs) ;



   % the new decision rule index is the capital index for which

   % vs is maximized



   drul(ikp,iap)=find(vs == vnew(ikp,iap)) ;

  end ;

 end ;



 % compute Tauchen's convergence criterion



 maxn=max(max(abs(vnew-vfn))') ;

 mind=min(min(vfn)') ;

 ccrit=maxn/mind ;



 iter=iter+1 ;

 [iter maxn mind ccrit]

 vfn=vnew ;



end ;



% compute the capitals that go with the decision rule index

drulk=reshape(kgrid(drul,1),nkp,nap) ;

% compute the capitals that correspond to the log-lin decision rule

drulkl=lkss+h1*(kgrid(:,ones(nap,1))-lkss)+h2*(agrid(:,ones(nkp,1))') ;



% print out some stuff





surf(agrid,kgrid,drulk)

title('Tauchen Solution')

xlabel('a(t)') ;

ylabel('k(t)') ;

zlabel('k(t+1)') ;



pause



surf(agrid,kgrid,drulkl)

title('Linearized Solution')

xlabel('a(t)') ;

ylabel('k(t)') ;

zlabel('k(t+1)') ;



pause



surf(agrid,kgrid,drulk-drulkl)

title('Difference')

xlabel('a(t)') ;

ylabel('k(t)') ;

zlabel(' ') ;



