function optionbinomial(S0,strike,sigma,T,N,type_op)
%Binomial model to price a stock option (European, American, or Bermudan)
%S0: current stock price, strike: option strike;
%sigma: volatility
%T: option expiration (in years)
%N: number of time steps
%type_op = 'call' or 'put'
type_ae = 'european';
type_md = 'sym_rw';
delta=0.0;               %dividend rate
dt=T/N;                  %time step
rf=0.01;                 %risk-free rate

% initialize stock price and option price arrays:
S=zeros(N,N+1);
V=zeros(N,N+1);

switch type_md
    case 'asym_rw'
        u=exp(sigma*sqrt(dt));   %up factor
        d=1.0/u;                 %down factor
        pu=0.5*(1+(rf-0.5*sigma*sigma)/sigma*sqrt(dt));
        pd=1-pu;
    case 'sym_rw'
        u = exp((rf-0.5*sigma*sigma)*dt+sigma*sqrt(dt));
        d = exp((rf-0.5*sigma*sigma)*dt-sigma*sqrt(dt));
        pu = 0.5;
        pd = 0.5;
end

disc = 1.0/(1.0+rf*dt); %exp(-rf*dt);      %discount factor in one time step

div_reduct = (1-delta*dt);

timetable=[10 20 30 40 50 60];

%set up the binomial tree for S:
S(1,1)=S0*d*div_reduct;
S(1,2)=S0*u*div_reduct;

for i=1:N-1
    S(i+1,1:i+1)=S(i,1:i+1)*d*div_reduct;
    S(i+1,i+2)=S(i,i+1)*u*div_reduct;
end

%generate the payoff function:

for j=1:N+1
    V(N,j) = payoff(S(N,j),strike,type_op);
end
%backward induction to obtain the option price:
switch type_ae
    case 'european'
    for i=N:-1:2
        V(i-1,1:i)=(pd*V(i,1:i)+pu*V(i,2:i+1))*disc;
    end
    case 'american'
    for i=N:-1:2
    for j=1:i
       V(i-1,j)=max(payoff(S(i-1,j),strike,type_op),(pd*V(i,j)+pu*V(i,j+1))*disc);
    end
    end
    case 'bermudan'
    for i=N:-1:2
        if any( (i-1)==timetable )
        for j=1:i
            V(i-1,j)=max(payoff(S(i-1,j),strike,type_op),(pd*V(i,j)+pu*V(i,j+1))*disc);
        end
        else
            V(i-1,1:i)=(pd*V(i,1:i)+pu*V(i,2:i+1))*disc;
        end
    end 
    otherwise
    fprintf('option style not recognized!\n');
end

switch type_ae
    case {'european','bermudan'}
        V0 = (pd*V(1,1)+pu*V(1,2))*disc;
    case 'american'
        V0 = max(payoff(S0,strike,type_op),(pd*V(1,1)+pu*V(1,2))*disc);
end

%generate the Black-Scholes price:
d1=(log(S0/strike)+(rf+0.5*sigma*sigma)*T)/(sigma*sqrt(T));
d2=d1-sigma*sqrt(T);
BS_put=-S0*0.5*erfc(d1/sqrt(2))+strike*exp(-rf*T)*0.5*erfc(d2/sqrt(2));
BS_call=S0*0.5*erfc(-d1/sqrt(2))-strike*exp(-rf*T)*0.5*erfc(-d2/sqrt(2));

fprintf('strike is %6.2f, vol is %4.2f%%, N = %d\n',strike,100*sigma,N);

switch type_op
    case 'call'
    fprintf('call price is %12.8f, BS price %10.6f, delta is %12.6f\n relative error %12.8f\n',...
             V0,BS_call,0.5*erfc(d1/sqrt(2)),(V0-BS_call)/BS_call);
    case 'put'
    fprintf('put price is %12.8f, BS price %10.6f, delta is %12.6f\n relative error %12.8f\n',...
             V0,BS_put,-0.5*erfc(-d1/sqrt(2)),(V0-BS_put)/BS_put);
    otherwise
    fprintf('option type not recognized!\n');
end
