% This program computes pricing of a European security derivative in a binomial
% pricing model.
%
% Radu Balan, January 28, 2008

% Input data
r = 0.25; % Interest rate
u = 2; % Up factor
d = 0.5; % Down factor
S0 = 4; % Stock initial share price
N = 3; % Number of periods
%VN = [27;3;3;0;3;0;0;0]; % VN: European call with strike price 5
VN = [0;8;0;6;0;2;2;3.5]; % VN: Lookback

% Compute risk-neutral probabilities
pt = (1+r-d)/(u-d);
qt = (u-1-r)/(u-d);

startind = zeros(N,1);
% Initialize price vector with VN
L = 2^N;
V = VN;
startind(N)=1; % Note we go backwards from end to the beginning
currentind = 1+L;
for k=N:-1:2
    % Allocate memory
    L = 2^(k-1);
    V = [V;zeros(L,1)];
    startind(k-1)=currentind;
    currentind = currentind + L;
    % Populate option price vector
    j1 = startind(k);
    j0 = startind(k-1);
    for k=1:L
        V(j0) = (pt*V(j1)+qt*V(j1+1))/(1+r);
        j0 = j0+1;
        j1 = j1+2;
    end
end
% Last value
V = [V;0];
lastind = currentind;
V(currentind) = (pt*V(startind(1)) + qt*V(startind(1)+1))/(1+r);
V0 = V(lastind);
D0 = (V(startind(1)) - V(startind(1)+1))/((u-d)*S0); % Initial short position

% Print prices
fprintf(1,' Option is sold at a premium of: %f at time zero\n',V0);
for k=1:N
    fprintf(1,'At time %d, option price and discounted option price are:\n',k)
    L = 2^k;
    for w=0:L-1
        fprintf(1,' %f  %f\n',V(startind(k)+w),V(startind(k)+w)/((1+r)^k));
    end
    fprintf(1,'\n\n');
end

% Prin portofolio evolution for a particular trajectory
% Market evolution
dynamics = 'HTH';
S = S0; % Current share price
X = V0; % Initial portofolio worth
D = D0; % Initial short position
M = X - D*S; % Initial position in money market
fprintf(1,'\nAt start time a portofolio worth %f composed of (%f shares, %f in money market)\n',X,D,M);
offset = 0;
for k=1:N-1
    switch dynamics(k),
        case 'H',
            S = S*u;
            offset = 2*offset;
        case 'T',
            S = S*d;
            offset = 2*(offset+1);
    end
    X = D*S + (1+r)*M;
    D = (V(startind(k+1)+offset) - V(startind(k+1)+offset+1))/((u-d)*S);
    M = X - D*S;
    fprintf(1,'At time %d a portofolio worth %f composed of (%f shares, %f in money market)\n',k,X,D,M);
%    fprintf(1,'   Control: Stock at %f per share, offset=%d\n\n',S,offset);
end
switch dynamics(N)
    case 'H',
        S = S*u;
    case 'T',
        S = S*d;
end
% Portofolio value
X = D*S + (1+r)*M;
fprintf(1,'End value of the portofolio: X=%f\n',X);