% 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