% Program that solves Ch 3 Exercise 2 (BCR part)

% 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

'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
% simulate the model to get std.dev. of log of K

nsim=1001 ;
a=normar1(0,rho,sigma,nsim) ;
k=zeros(nsim,1) ;

ccrit=1 ;
sigm=zeros(2,2) ;
sige=[ 0.0 0.0 ; 0.0 (sigma*sigma) ]
mmat=[ h1 h2 ; 0.0 rho ]
gmat=eye(2) ;
while (ccrit > 1e-12) ;
 badd=gmat*sige*(gmat') ;
 ccrit=max(max(abs(badd))') ;
 sigm=sigm+badd ;
 gmat=mmat*gmat ;
end ;
sdk=sqrt(sigm(1,1))

% set up the grid for k

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

% initialize the decision rule by obtaining the decision rule closest
% to the log-linear case

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 decision rule

drul=zeros(nkp,nap) ;
drulk=zeros(nkp,nap) ;
for ikp=1:nkp ;
 for iap=1:nap ;

  % find the k in the grid that's closest to the log-linear rule

  gap=abs(kgrid-h1*kgrid(ikp,1)-h2*agrid(iap,1)) ;
  drul(ikp,iap)=find(gap == min(gap)) ;
  drulk(ikp,iap)=kgrid(drul(ikp,iap),1) ;
 end ;
end ;

lkss=log(kss) ;    % log of nssv of K for A=1 case
kgrid=kgrid+lkss ; % add the mean to the k grid
drulk=drulk+lkss ; % add the mean to the initial decision rule
drulkl=drulk ;     % store this initial guess permanently

bigix=cumsum(ones(nkp,nap)) ; % creates a matrix whose ji element is j

% now iterate on the decision rule to a fixed point

iter=0 ;
ccrit=1 ;
while ((iter <= 10000)&(ccrit > 0)) ;

 dnew=zeros(nkp,nap) ;  % initialize the new decision rule indexes
 dnewk=zeros(nkp,nap) ; % initialize the new decision rule k values

 for ikp=1:nkp ;
  for iap=1:nap ;

   % evaluate the piece on the left-hand side

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

   lhs=ones(nkp,1)./ccurp ;

   % evaluate the piece on the right-hand side

   rfut=thet*exp((agrid(:,ones(nkp,1))')+(thet-1).* ...
        kgrid(:,ones(nap,1)))+(1-delt) ;
   cfut=exp((agrid(:,ones(nkp,1))')+thet.*kgrid(:,ones(nap,1)))+ ...
        (1-delt)*exp(kgrid(:,ones(nap,1)))-exp(drulk) ;
   rhs=bet.*((rfut./cfut)*(pim(iap,:)')) ;

   gap=abs(lhs-rhs) ; % evaluate the difference

   dnew(ikp,iap)=find(gap == min(gap)) ;   % find where they are close
   dnewk(ikp,iap)=kgrid(dnew(ikp,iap),1) ;
  end ;
 end ;
 ccrit=sum(sum(dnew ~= drul)') ; % counts number of points at which
                                 % the decision rule changed

 % the next statement finds out if there are any attracting states
 % for capital, i.e. any points at which k(i) is optimal at k(i)
 % and all a(j) (it counts the number of such points)

 trapk=  sum(  (sum((drul == bigix)')') == nap) ;

 iter=iter+1 ;
 [iter ccrit trapk]
 drul=dnew ;   % update the rules
 drulk=dnewk ;

end ;

diary drules.out ;
[drulk(:,1) drulkl(:,1)]
[drulk(:,nap) drulkl(:,nap)]
diary off ;
plot(kgrid,[drulk drulkl])
pause
plot(kgrid,drulk-drulkl)
