%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
%       ayse.m: replicates Ayse Imrohoroglu's JEDC paper
%
%   George Hall
%   Yale University
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
clear
! rm ayse.out
diary ayse.out
disp('AYSE MODEL');
disp('');
t1 = clock;
%
%  set parameter values
%
sigma  = 1.50;                % risk aversion              
beta   = 0.995;               % subjective discount factor 
prob = [ .9565 .0435; .5 .5]; % prob(i,j) = probability(s(t+1)=sj|s(t)=si)
tau    = 0.0125;              % growth rate of money supply
theta  = 0.25;                % non-interest income if unemploy is theta*wage
wage   = 1.00;                % non-interest income if employed
Astart = 2.267;               % initial value for aggregate level of assets
g      = 1.00;                % relaxation parameter
R      = 1/(1+tau);           % gross rate of interest
w      = 0.0;                 % compensating variation           
%
%   form asset grid
%   
maxast =  8.0;                     % maximum value of asset grid   
incast = .026667;                  % size of asset grid increments
nasset = round(maxast/incast+1);   % number of grid points
agrid = [ 0:.026667: 300*incast ]';        % asset grid
%
%   loop to find fixed point for agregate level of assets, A. 
%
liter   = 1;
maxiter = 40;
toler   = 1e-4;
metric  = 100;
A = Astart;
disp('ITERATING ON A');
disp('');
disp('    liter     metric     meanA     Aold');
while  (metric > toler) & (liter <= maxiter);
  %
  %  tabulate the utility function such that for zero or negative
  %  consumption utility remains a large negative number so that
  %  such values will never be chosen as utility maximizing      
  %
  util1=-inf*ones(nasset,nasset);  % utility when employed     
  util2=-inf*ones(nasset,nasset);  % utility when unemployed   
  for i=1:nasset;
    asset=(i-1)*incast;
    for j=1:nasset; 
      assetp = (j-1)*incast;
      cons = w + wage + R*asset + A*R*tau - assetp;
      if cons > 0; 
	util1(j,i)=((cons)^(1-sigma)-1)/(1-sigma);
      end;
    end;
    for j=1:nasset
      assetp = (j-1)*incast;
      cons = w + theta*wage + R*asset + A*R*tau - assetp;
      if cons > 0;
	util2(j,i)=((cons)^(1-sigma)-1)/(1-sigma);
      end;
    end;
  end;
  %
  %  initialize some variables
  %
  v       = zeros(nasset,2);
  decis   = zeros(nasset,2);
  test1    = 10;
  test2    = 10;
  [rs,cs] = size(util1);
  %
  %  iterate on Bellman's equation and get the decision 
  %  rules and the value function at the optimum         
  %
  while (test1 ~= 0) | (test2 >= 1e-2);
    for i=1:cs;
      r1(:,i)=util1(:,i)+beta*(prob(1,1)*v(:,1)+ prob(1,2)*v(:,2));
      r2(:,i)=util2(:,i)+beta*(prob(2,1)*v(:,1)+ prob(2,2)*v(:,2));
    end;
    
    [tv1,tdecis1]=max(r1);
    [tv2,tdecis2]=max(r2);
    tdecis=[tdecis1' tdecis2'];
    tv=[tv1' tv2'];
    
    test1=max(any(tdecis-decis));
    test2=max(max(abs(tv-v))'); 
    v=tv;
    decis=tdecis;
  end;
  decis=(decis-1)*incast;
  %
  %   form transition matrix
  %   trans is the transition matrix from state at t (row)
  %   to the state at t+1 (column) 
  % 
  g2=zeros(cs);
  g1=zeros(cs);
  for i=1:cs
    g1(i,tdecis1(i))=1;
    g2(i,tdecis2(i))=1;
  end
  trans=[ prob(1,1)*g1 prob(1,2)*g1; prob(2,1)*g2 prob(2,2)*g2];
  trans= trans';
  strans = sparse(trans);
  probst = (1/(2*nasset))*ones(2*nasset,1);
  test = 1;
  while test > 10^(-8);
    probst1 = strans*probst;
    test=max(abs(probst1-probst));
    probst = probst1;
  end;
  %
  %   vectorize the decision rule to be conformable with probst
  %   calculate new aggregate level of assets   meanA
  %
  aa=decis(:);
  meanA=probst'*aa;
  %
  %  calculate measure over (k,s) pairs
  %  lambda has same dimensions as decis
  %
  lambda=zeros(cs,2);
  lambda(:)=probst;
  %
  %   calculate stationary distribution of a 
  %
  probk=sum(lambda');     %  stationary distribution of `captal' 
  probk=probk';
  %
  %   form metric and update A
  %
  Aold = A;
  Anew = g*meanA + (1-g)*Aold;
  metric = abs((Aold-meanA)/Aold);
  A = Anew;
  disp([ liter metric meanA Aold ]);
  liter = liter+1;
end;
%timer = etime(clock,t1);
%disp('Finding the Fixed Point took:');
%disp([timer]);
%disp('seconds');
%
%   Calculate Expected Utility
%
meanA2  = probst'*(decis(:).^2);
varA    = meanA2 - meanA^2;
stdA    = sqrt(varA);
grid    = linspace(0,maxast,nasset);
congood = w + wage + R*grid' + A*R*tau - grid(tdecis(:,1))';
conbad  = w + theta*wage + R*grid' + A*R*tau - grid(tdecis(:,2))';
consum  = [congood conbad ];
cons2   = [congood.^2 conbad.^2];
meancon = sum(diag(lambda'*consum));
meancon2= sum(diag(lambda'*cons2));
varcon  = ( meancon2 - meancon^2 );
stdcon  = sqrt(varcon);
UCEU    = sum(diag(lambda'*v));
UTILITY = (consum.^(1-sigma)-1)./(1-sigma);
UCEU2   = sum(diag(lambda'*UTILITY));
%
%   print out results
%
disp('PARAMETER VALUES');
disp('');
disp('    sigma      beta        R       tau     theta    w'); 
disp([ sigma beta R tau theta w]);
disp(''); 
disp('ASSET GRID');
disp('');
disp('    maxast     incast');
disp([ maxast incast ]);
disp('');
disp('EQUILIBRIUM RESULTS ');
disp('');
disp('      A       stdA    UCEU    UCEU2      MEANCON     STDCON');
disp([ Aold stdA UCEU UCEU2 meancon stdcon]);

%
%   plot invariant distribution of money holdings 
%

figure(1)
plot(agrid,probk);
xlabel('money holdings')
ylabel('fraction of agents');
print mondist.ps

exit