function [fp dyn] = fnct(eta, f, r, rho, a, a_, g, delta_, sigma)
% [fp dyn] = fnct(eta, f, r, rho, a, a_, g, delta_, sigma)
% takes scalar eta, a 4x1 vector f = [theta, theta', q, q'], and parameters
% r, rho, a, sigma, g, delta_, and computes fp, the derivative of f with
% respect to eta, as well as dyn = [psi, sigma_eta*eta, sigma_q, mu_eta*eta, mu_q]
%
% assumes that for both experts and households, the production sets consist
% of single points (a, g) and (a, - delta_) respectively 
%
% written by Yuliy Sannikov

    

% search for psi between eta (lower bound) and min(f(3)/f(4) + eta, 1) (upper bound) 
    
psi_L = eta; psi_R = min(f(3)/f(4) + eta, 1);  % alpha
for n = 1:30,
            psi = (psi_L + psi_R)/2;
            amplification = 1 - f(4)/f(3)*(psi - eta);
            
            % VOLATILITY COMPUTATION
            sigma_eta_eta = sigma*(psi - eta)/amplification;  % sigma_eta *times* eta
            sigma_q = sigma_eta_eta*f(4)/f(3); 
            sigma_theta = sigma_eta_eta*f(2)/f(1);
            risk_premium = - sigma_theta*(sigma + sigma_q);
            
            household_premium = (a_ - a)/f(3) - delta_ - g + risk_premium;
                   
            if household_premium > 0, % households want to hold more
                psi_R = psi;
            else
                psi_L = psi;
            end
end

mu_q = r - a/f(3) - g - sigma*sigma_q + risk_premium;
mu_eta_eta = - (psi - eta)*(sigma + sigma_q)*(sigma + sigma_q + sigma_theta) + eta*a + eta*(1 - psi)*(delta_ + g);
qpp = 2*(mu_q*f(3) - f(4)*mu_eta_eta)/sigma_eta_eta^2; 
thetapp = 2*((rho - r)*f(1) - f(2)*mu_eta_eta)/sigma_eta_eta^2; 

fp = [f(2); thetapp; f(4); qpp];

dyn = [psi, sigma_eta_eta, sigma_q, mu_eta_eta, mu_q];





