%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%% Delta-Gamma Hedging using time and moved based rebalancing times 
%%% incorporating transaction costs.
%%% 
%%% Copyright (C) 2009  Sebastian Jaimungal
%%%
%%%
%%% This program is free software: you can redistribute it and/or modify
%%% it under the terms of the GNU General Public License as published by
%%% the Free Software Foundation, either version 3 of the License, or
%%% (at your option) any later version.
%%%
%%% This program is distributed in the hope that it will be useful,
%%% but WITHOUT ANY WARRANTY; without even the implied warranty of
%%% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
%%% GNU General Public License for more details.
%%% 
%%% You should have received a copy of the GNU General Public License
%%% along with this program.  If not, see <http://www.gnu.org/licenses/>.
%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function DGHedge

    Nsims = 50;
    
    % turn to true to plot sample paths
    % turn to false to only show histograms -- reccommended when Nsims is
    % large
    DoSamplePathPlot = true;
    
    % size of moves in delta before rebalancing
    DDelta = 0.05;

    % transaction costs
    transS =  0.001;
    transOpt =  0.005;

    % maturity and number of steps
    T = 1;
    Nsteps = 252;
    
    % model parameters
    mu = 0.2;
    volP = 0.2;
    volI = 0.2;
    S0 = 100;
    K = 100;
    rf = 0.01;
    
    
    % grab sample paths
    Spot = SimPath(Nsims, Nsteps, T, S0, mu, volP);
    
    dt = T/Nsteps;
    
    % matrix of time to maturities
    Terms = T - dt*(cumsum(ones(Nsims, Nsteps+1),2) - 1)+1e-10;
    
    % compute Deltas and Gammas with given sample path
    DeltaTimeBased = PutDelta(Spot, K, Terms, rf, volI);
    DeltaMoveBased = DeltaTimeBased;
    
    GammaTimeBased = PutGamma(Spot, K, Terms, rf, volI);
    
    % compute portfolio value through time for each sample path
    PutPrices = PutPrice(Spot, K, Terms, rf, volI);
    
    % compute hedging portfolio
    FPrice = PutPrice(Spot, K, Terms+0.25, rf, volI);
    FDelta = PutDelta(Spot, K, Terms+0.25, rf, volI);
    FGamma = PutGamma(Spot, K, Terms+0.25, rf, volI);
    
    NumFTimeBased =  GammaTimeBased ./ FGamma;
    NumFMoveBased = NumFTimeBased;
    
    NumSTimeBased = DeltaTimeBased - NumFTimeBased .* FDelta;
    NumSMoveBased = NumSTimeBased;
   
        
    % store bank account amounts using time-based and move-based hedging
    BankMoveBased = zeros(Nsims, Nsteps+1);
    BankMoveBased(:,1) = PutPrices(1,1) - Spot(1,1) * NumSMoveBased(1,1) - NumFMoveBased(1,1) * FPrice(1,1) ...
                         - abs(NumSMoveBased(1,1)) * Spot(1,1) * transS ...
                         - abs(NumFMoveBased(1,1)) * FPrice(1,1) * transOpt;
    
    BankTimeBased = BankMoveBased;
    
    % store levels surrounding current delta hedge position
    DeltaBd = zeros(Nsims, 2);
    
    % stores whether a trade is made in move-based approach
    KeepDelta = zeros(Nsims,Nsteps);
    
    for i = 2 : Nsteps+1

        % compute bounding delta-levels
        DeltaBd(:,1) = ( NumSMoveBased(:,i-1) + NumFMoveBased(:,i-1) .* FDelta(:,i-1) ) - DDelta;
        DeltaBd(:,2) = DeltaBd(:,1) + 2 * DDelta;

        
        % stores whether or not the delta has moved enough
        RepTimeDelta = NumSTimeBased(:,i) + NumFTimeBased(:,i) .* FDelta(:,i);
        
        KeepDelta(:,i-1) =     ( RepTimeDelta > DeltaBd(:,1) ) ...
                            .* ( RepTimeDelta < DeltaBd(:,2));
        
        % update positions if moved enough, otherwise keep old positions
        NumSMoveBased(:,i) =   NumSMoveBased(:,i-1) .* KeepDelta(:,i-1) ...
                              + NumSTimeBased(:,i) .* (1-KeepDelta(:,i-1));
                          
        NumFMoveBased(:,i) =   NumFMoveBased(:,i-1) .* KeepDelta(:,i-1) ...
                              + NumFTimeBased(:,i) .* (1-KeepDelta(:,i-1));                          
        
        DeltaMoveBased(:,i) = NumSMoveBased(:,i) + NumFMoveBased(:,i) .* FDelta(:,i);
                          
        % update bank account
        BankMoveBased(:,i) =  BankMoveBased(:,i-1)*exp(rf*dt) ...
                            + ( NumSMoveBased(:,i-1) - NumSMoveBased(:,i) ) .* Spot(:,i) ...
                            - abs( NumSMoveBased(:,i-1) - NumSMoveBased(:,i) ) .* Spot(:,i) * transS ...
                            + ( NumFMoveBased(:,i-1) - NumFMoveBased(:,i) ) .* FPrice(:,i) ...
                            - abs( NumFMoveBased(:,i-1) - NumFMoveBased(:,i) ) .* FPrice(:,i) * transOpt;
                        
        BankTimeBased(:,i) =  BankTimeBased(:,i-1)*exp(rf*dt) ...
                            + ( NumSTimeBased(:,i-1) - NumSTimeBased(:,i) ) .* Spot(:,i) ...
                            - abs( NumSTimeBased(:,i-1) - NumSTimeBased(:,i) ) .* Spot(:,i) * transS ...
                            + ( NumFTimeBased(:,i-1) - NumFTimeBased(:,i) ) .* FPrice(:,i) ...
                            - abs( NumFTimeBased(:,i-1) - NumFTimeBased(:,i) ) .* FPrice(:,i) * transOpt; 
        
    end
    
    % compute total pnl through the entire sample path
    pnlTimeBased = Spot.* NumSTimeBased + FPrice .* NumFTimeBased + BankTimeBased - PutPrices;
    pnlMoveBased = Spot.* NumSMoveBased+ FPrice .* NumFMoveBased +  BankMoveBased - PutPrices;
    

    % plot results....
    
    close all;
    
    if DoSamplePathPlot
        
        % Time Based...
        figure

        subplot(4,1,1);
        plot(Spot');
        %xlabel('time', 'Fontsize',12)
        ylabel('Spot Price', 'Fontsize',12)
        title('Time Based -- Spot price paths', 'Fontsize',12)
        axis([0 Nsteps+1 S0*0.5 S0*1.5])

        subplot(4,1,2);
        plot(pnlTimeBased' )
        %xlabel('time', 'Fontsize',12)
        ylabel('Value', 'Fontsize',12)
        title('Time Based -- Hedged Portfolio Value', 'Fontsize',12)
        axis([0 Nsteps+1 -2 2])


        subplot(4,1,3);
        plot(DeltaTimeBased')
        xlabel('time', 'Fontsize',12)
        ylabel('Delta', 'Fontsize',12)
        title('Time Based -- Delta Position', 'Fontsize',12)
        axis([0 Nsteps+1 -1 0]);
        
        subplot(4,1,4);
        plot((NumFTimeBased.*GammaTimeBased)')
        xlabel('time', 'Fontsize',12)
        ylabel('Gamma', 'Fontsize',12)
        title('Time Based -- Gamma Position', 'Fontsize',12)        
        axis([0 Nsteps+1 0 0.5]);
        
        % Move Based
        figure

        subplot(4,1,1);
        plot(Spot');
        %xlabel('time', 'Fontsize',12)
        ylabel('Spot Price', 'Fontsize',12)
        title('Move Based -- Spot price paths', 'Fontsize',12)
        axis([0 Nsteps+1 S0*0.5 S0*1.5])

        subplot(4,1,2);
        plot(pnlMoveBased' )
        %xlabel('time', 'Fontsize',12)
        ylabel('Value', 'Fontsize',12)
        title('Move Based -- Hedged Portfolio Value', 'Fontsize',12)
        axis([0 Nsteps+1 -2 2])

        subplot(4,1,3);
        plot(DeltaMoveBased')
        xlabel('time', 'Fontsize',12)
        ylabel('Delta', 'Fontsize',12)
        title('Move Based -- Delta Position', 'Fontsize',12)
        axis([0 Nsteps+1 -1 0]);
        
        subplot(4,1,4);
        plot((NumFMoveBased.*GammaTimeBased)')
        xlabel('time', 'Fontsize',12)
        ylabel('Gamma', 'Fontsize',12)
        title('Move Based -- Gamma Position', 'Fontsize',12)        
        axis([0 Nsteps+1 0 0.5]);        
      
        
        % Positions comparison
        figure
        
        subplot(4,1,1);
        plot(NumSTimeBased')
        xlabel('time', 'Fontsize',12)
        ylabel('Num S', 'Fontsize',12)
        title('Time Based -- S position', 'Fontsize',12)
        
        subplot(4,1,2);
        plot(NumSMoveBased')
        xlabel('time', 'Fontsize',12)
        ylabel('Num S', 'Fontsize',12)
        title('Move Based -- S position', 'Fontsize',12)
        
        subplot(4,1,3);
        plot(NumFTimeBased')
        xlabel('time', 'Fontsize',12)
        ylabel('Num F', 'Fontsize',12)
        title('Time Based -- F position', 'Fontsize',12)
        
        subplot(4,1,4);
        plot(NumFMoveBased')
        xlabel('time', 'Fontsize',12)
        ylabel('Num F', 'Fontsize',12)
        title('Move Based -- F position', 'Fontsize',12)        
        
    end
        
    % historgram of terminal PnL
    bins = [-5:0.25:5];
    figure
    [ freq bins ] = hist(pnlTimeBased(:,Nsteps+1), bins);
    hist(pnlTimeBased(:,Nsteps+1), bins);
    axis([bins(1) bins(end) 0 max(freq)*1.1 ]);
    title('Time Based P n L');
    
    figure
    [ freq bins ] = hist(pnlMoveBased(:,Nsteps+1), bins);
    hist(pnlMoveBased(:,Nsteps+1), bins);
    axis([bins(1) bins(end) 0 max(freq)*1.1 ]);
    title('Move Based P n L');
    
    % compare replication with payoff
    terminalReplTimeBased = Spot(:,end) .* NumSTimeBased(:,end) + FPrice(:,end) .* NumFTimeBased(:,end) + BankTimeBased(:,end);
    terminalReplMoveBased = Spot(:,end) .* NumSMoveBased(:,end) + FPrice(:,end) .* NumFMoveBased(:,end) + BankMoveBased(:,end);
    
    figure
    scatter(Spot(:,Nsteps+1), terminalReplTimeBased)
    hold
    plot( [0:1:250], max( K - [0:1:250], 0), '-r','LineWidth',3);
    hold
    title('Time Based Replication');
    axis([ 0 250 -25 125]);    
    
    figure
    scatter(Spot(:,Nsteps+1), terminalReplMoveBased)
    hold
    plot( [0:1:250], max( K - [0:1:250], 0), '-r','LineWidth',3);
    hold
    title('Move Based Replication');
    axis([ 0 250 -25 125]);    
    
    % plot histogram of number of trades using moved-based trading
    figure
    Ntrades = sum(1-KeepDelta,2);
    hist(Ntrades, [0:1:255]);
    title('Histogram of Number of Move-Based Trades');
    
end

function price = CallPrice(S, K, T, rf, sigma)

    dp = (log(S/K) + (rf + 1/2*sigma.^2).*T) ./ (sigma.*sqrt(T));
    dm = dp - sigma.*sqrt(T);
    
    price = S .* normcdf(dp, 0, 1) - K .* exp(-rf*T) .* normcdf(dm, 0, 1);
    
end

function delta = CallDelta(S, K, T, rf, sigma)

    dp = ( log(S ./ K) + (rf + 1/2 * sigma.^2).* T ) ./ (sigma.*sqrt(T));
    delta = normcdf(dp, 0, 1);
    
end

function gamma = CallGamma(S, K, T, rf, sigma)

    dp = ( log(S./K) + (rf + 1/2*sigma.^2).*T) ./ (sigma.*sqrt(T));
    gamma = normpdf(dp, 0, 1) ./ (sigma .* S .* sqrt(T) );
    
end

function vega = CallVega(S, K, T, rf, sigma)

    dp = ( log(S./K) + (rf + 1/2*sigma.^2).*T) ./ (sigma.*sqrt(T));
    
    vega = S .* normpdf(dp, 0, 1) .* sqrt(T); 
    
end


function price = PutPrice(S, K, T, rf, sigma)

    dp = ( log(S./K) + (rf + 1/2*sigma.^2).*T) ./ (sigma.*sqrt(T));
    dm = dp - sigma.*sqrt(T);
    
    price =  - S .* normcdf(-dp, 0, 1) + exp(-rf*T) .* K .* normcdf(-dm, 0, 1);
    
end

function delta = PutDelta(S, K, T, rf, sigma)

    dp = (log(S./K) + (rf + 1/2*sigma.^2).*T) ./ (sigma.*sqrt(T));
    delta = normcdf(dp, 0, 1) - 1;
    
end

function gamma = PutGamma(S, K, T, rf, sigma)

    dp = (log(S./K) + (rf + 1/2*sigma.^2).*T) ./ (sigma.*sqrt(T));
    gamma = normpdf(dp, 0, 1) ./ (sigma .* S .* sqrt(T) );
    
end

function vega = PutVega(S, K, T, rf, sigma)

    dp = (log(S./K) + (rf + 1/2*sigma.^2).*T) ./ (sigma.*sqrt(T));
    
    vega = S .* normpdf(dp, 0, 1) .* sqrt(T); 
    
end