Pricing American Put Options via Binomial Tree in Matlab

Question

I currently am completing a Computational Finance Assignment, and am trying to figure out how to alter this Matlab code which prices a European put or call option, in order to price an American Put Option. I honestly thought it would be as simple as placing a max() in the backwards recursion step. I don't want you to just provide the altered code, as I'd rather learn, but I have been thinking about this for a while and am at crossroads. I have left my altered code in thus far in the hope that you could point me in the right direction.

function price  = tree_slow(S0, K, T, r, sigma, opttype, Nsteps)
%   
% S0 - current stock price
% K - strike
% T - expiry time
% r - interest rate
% sigma - volatility
% opttype - 0 for a call, otherwise a put
% Nsteps - number of timesteps

%Output
%   price : option price

%Practical 1: compute the timestep size (Delta t) and tree parameters
delt = T/Nsteps;
u = exp(sigma * sqrt(delt) );
d = 1./u;
a = exp( r*delt );
p = (a - d)/(u - d);

%vector of payoff and option price in the tree
W = zeros(Nsteps+1,1);

%Practical 1: compute the S value at time T and store it in W
for j=0:Nsteps
    W(j+1,1) = S0*u^(j)*d^(Nsteps -j);
end

%Practical 1: compute the payoff
if(opttype == 0)
    W = max(W-K,0); 
else
    W = max(K-W,0);
end

%Practical 2: fill in the backward recursion
for n=Nsteps-1:-1:0%timeloop

    %loop over all possible S levels at time t_n 
    for j=0:n 
        %instruction: complete the expectation formula
        W(j+1,1) = max(K-W(j+2,1),exp(-r*delt)*( p*W(j+2,1) + (1-p)*W(j+1,1) ));
    end

end
%instruction: fill in with the right index
price = W(1);

Answer 1

It is as simple as just taking the max(). The problem is that you took the wrong one.

You must consider the max between the intrinsic value of the option on the one hand and its discounted continuation value (which is an expectation in the risk-neutral world) on the other.

In your final loop, you should therefore replace the line

W(j+1,1) = max(K-W(j+2,1),exp(-r*delt)*( p*W(j+2,1) + (1-p)*W(j+1,1) ));

with

W(j+1,1) = max(phi, exp(-r*delt)*( p*W(j+2,1) + (1-p)*W(j+1,1) ));

where for the tree node j at the time iteration n

w = (+1) * (opttype == 0) + (-1) * (opttype ~= 0)    % w = (opttype == 0) ? +1 : -1
S = S0*u^(j)*d^(n-j)
phi = w*(S-K)