% This file computes the equilibrium in the setting of Brunnermeier and
% Sannikov, under assumptions that (1) experts can only issue equity and 
% (2) the production sets of both experts and households consist
% of single points (a, g) and (a_, -delta_) respectively 
%
% written by Yuliy Sannikov

tic

color = 'g';
global qmax
a = 1; a_ = 0.5; rho = 0.06; r = 0.05; sigma = 0.01; g = 0.04; delta_ = 0.05;
q_ = a_/(r + delta_);
qmax = a/(r - g);


etaspan = [0 1];
F0 = [1 -1e+10 q_ 0]';  % the last coordinate is the derivative q'(0)
                        % found in the following loop

options = odeset('events','evntfcn');  % this function terminates integration
         % of the ode when q exceeds qmax, or f'(eta) or q'(eta) reaches 0 
                                    
odefun = @(eta,f) fnct(eta, f, r, rho, a, a_, g, delta_, sigma);

QL = 0; QR = 1e+15;
for iter = 1:50,
    F0(4) = (QL + QR)/2;  % this is q'(0)
    % solves the system of ODE's starting with boundary conditions 
    % theta'(0) = #large negative number# and q(0) = q_; searching 
    % for q'(0) such that q'(eta*) = theta'(eta*) at the same level of eta*
    % note that if theta(eta) satisfies the system of ODE's then so does
    % const*theta(eta), for any positive constant const

    [etaout,fout,TE,YE,IE] = ode45(odefun,etaspan,F0,options);
    
    if IE == 3, % if q'(eta) has reached zero, we 
        QL = F0(4);  % increase q'(0)
    else
        QR = F0(4);  % else we reduce q'(0)
    end

end
    
N = length(etaout);
Dyn = zeros(N, 5);

% calculate the psi, as well as drifts and volatilities 
for n = 1:N,
    [fp Dyn(n,:)] = fnct(etaout(n), fout(n,:), r, rho, a, a_, g, delta_, sigma);
end

% normalize theta, to make sure that theta(eta*) = 1
normalization = fout(N,1);
fout(:,1:2) = fout(:,1:2)/normalization;

figure(2);

subplot(3,2,1); hold on
plot(etaout(1:N-1), fout(1:N-1,3), color);   
xlabel('\eta')
ylabel('q');

subplot(3,2,2); hold on
plot(etaout(1:N-1), fout(1:N-1,1), color);   
xlabel('\eta')
ylabel('\theta');

subplot(3,2,3); hold on
plot(etaout(1:N-1), Dyn(1:N-1, 4), color);  
ylabel('\eta \mu_\eta');


subplot(3,2,4); hold on
plot(etaout(1:N-1), Dyn(1:N-1, 2), color);   
ylabel('\eta \sigma_\eta');
    
        
subplot(3,2,5); hold on
plot(etaout(1:N-1), Dyn(1:N-1, 1), color);   
ylabel('\psi');
        
subplot(3,2,6); hold on
plot(etaout(1:N-1), Dyn(1:N-1,3), color);   
ylabel('\sigma^q');


toc