function [V, sigmaQ, muQ] = ybinom_eu(rc,S0,sigma,mu,T,H,N)
% [V0, sigmaQ, muQ] = ybinom_eu(rc,S0,sigma,mu,T,H,N)
% find European option price using binomial tree method using yearly data rc,sigma,mu
% INPUT:
%     rc:  yearly interest rate for continuous compounding
%     S0:  initial stock price
%  sigma:  yearly volatility of stock price
%     mu:  yearly drift of stock price
%      T:  time to maturity in years
%      H:  payoff function, e.g. @(s)max(5-s,0) for put option with strike 5
%      N:  number of time steps for binomial tree model (period = T/N)
% OUTPUT:
%     V0:  European option price at time t=0
% sigmaQ:  yearly volatility w.r.t. Q
%    muQ:  yearly drift w.r.t. Q
%
% example:
%   rc = 0.10;              % interest rate 10%
%   H = @(s) max(s-7,0);    % call option with strike K=7 
%   T = 1.5;                % maturity: 1.5 years from now
%   S0 = 6; sigma = 0.5; mu = 0.8;
%   
%   for N = 2.^(5:10)
%     [V0,sigmaQ,muQ] = ybinom_eu(rc,S0,sigma,mu,T,H,N);
%     disp([N,V0,sigmaQ,muQ])
%   end   

n = N/T;                             % number of periods per year
rho = exp(rc/n) - 1;                 % interest rate per period 1/n
alpha = sqrt(sigma^2/n+mu^2/n^2);
u = exp(alpha); d=1/u;

q = (1+rho-d)/(u-d);                 % risk-neutral probability
epsilon = q - .5;
muQ = 2*n*alpha*epsilon;
sigmaQ = sqrt( n*alpha^2*(1-4*epsilon^2) );

S = S0*d^N*(u/d).^(0:N);             % stock prices at time t_N
V = H(S);                            % option prices at time t_N given by payoff function H
for j=N-1:-1:0                       % find option prices at time t_j, j=N-1,...,0
  V = (q*V(2:end)+(1-q)*V(1:end-1))/(1+rho);    % risk-neutral formula for option price
end

