%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%   Bastian Groß
%   Universität Trier
%
%   Numerische Lösung der Black-Scholes-Gleichung mit dem
%   Finite-Differenzen-Verfahren
%
%   Nach der Variablentransformation tau = T - t sowie x = ln(K)
%   ergibt sich die Differentialgleichung:
%
%   u_tau = sigma^2/2 * u_xx + (r - sigma^2/2) * u_x - r * u   auf [0,T] x (-G, G)
%       Anfangsbedingung:   u(x,0) = max(exp(x) - B, 0)       für alle x aus (-G,G)
%       Randbedingung:      u(-G,tau) = 0                     für alle tau aus [0,T]
%                           u(G,tau) = exp(G) - B*exp(-tau*r) für alle tau aus [0,T]
%
%   Die numerische Lösung wird mittels der geschlossenen Lösungsformel mit
%   exakten Werten verglichen
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%


function FD_Black_Scholes(x_disk, tau_disk, theta)



% ---------------------------------------
% ------------- Konstanten --------------
% ---------------------------------------
    T =     0.5;    % Zeithorizont
    G =     5;      % legt das Intervall [-G,G] fest
    sigma = 0.3;    % Volatilität
    r =     0.03;   % risikoloser Zinssatz 
    B =     50;     % Strike-Preis (Ausübepreis)
    
    % Diskretisierungsparameter
    N = x_disk;     % Anzahl der Schritte in x-Richtung
    M = tau_disk;   % Anzahl der Schritte in tau-Richtung
    delta_x   = (2*G) / (N+1);
    delta_tau = T / (M-1);
    hilf_x    = [-G + delta_x : delta_x : G - delta_x]';
    hilf_tau  = [0 : delta_tau : T]';    
    
% --------------------------------------
% ---- Initialisierung der Matrizen ----
% --------------------------------------
    % Lösungsmatrix "U"
    U = zeros(N, M);    % Jede Spalte U(:,i) entspricht dem Lösungsvektor in einem Zeitschritt
   
    % Startvektor
    U(:,1) = max(exp(hilf_x)-B,0);
    
    % Steifigkeitsmatrix "A"
    alpha =  sigma^2   + delta_x^2 * r;                % Hauptdiagonale
    beta  = -sigma^2/2 - delta_x/2 * (r - sigma^2/2);  % obere Nebendiagonale
    gamma = -sigma^2/2 + delta_x/2 * (r - sigma^2/2);  % untere Nebendiagonale
        
    e = ones(N,1);     % Hilfsvektor
    A = spdiags([gamma*e alpha*e beta*e], -1:1, N,N);
    
   % Randvektor "f"
    f_i_1 = zeros(N,1);  % f_i_1 steht für f im Zeitschritt (i-1)
    % f(1) = (delta_tau / delta_x^2) * ( sigma^2/2 - delta_x/2 * (r - sigma^2/2)) *  0;      % = 0!
    f_i_1(N) = (delta_tau / delta_x^2) * ( sigma^2/2 + delta_x/2 * (r - sigma^2/2)) *  (exp(G)- B);
    
    
% --------------------------------------
%  Lösen der linearen Gleichungssystems 
% --------------------------------------    

    % Allgemeine theta-Methode:
    % theta = 1: implizitere Euler, theta = 0: expliziter Euler, theta = 0.5: Crank-Nicolson
    for i = 2 : M
       % Neuer Randvektor (da dieser zeitabhängig ist!)
        f_i = zeros(N,1);
        f_i(N) = (delta_tau / delta_x^2) * ( sigma^2/2 + delta_x/2 * (r - sigma^2/2)) *  (exp(G)- B*exp(-hilf_tau(i)*r));   %B*exp(-(i-1)*delta_tau*r));
       % Berechnen der neuen Lösung
        U(:,i) = (eye(N) + (theta) * (delta_tau / delta_x^2) * A) \  ( ( eye(N) - (1-theta) * (delta_tau / delta_x^2) * A) * U(:,i-1) + theta*f_i + (1-theta)*f_i_1 );
       % Speichern des Vektors f_i als f_i_1
        f_i_1 = f_i;
    end
   
    
% --------------------------------------
%  Berechnen der Lösung mittels der 
%  geschlossenen Lösungsformel "blsprice"
% --------------------------------------        
    for i = 1 : M
        U_exakt(:,i) = blsprice(exp(hilf_x), B, r, hilf_tau(i), sigma); 
    end
        
% --------------------------------------
%  Darstellung des Fehlers zwischen 
%  Finite-Differenzen und BS
% --------------------------------------       
    figure(1);
    surf( T - hilf_tau,  exp(hilf_x), U - U_exakt)    
   
    figure(2);
    surf( T - hilf_tau,  exp(hilf_x), U )
    
    figure(3);
    surf( T - hilf_tau,  exp(hilf_x), U_exakt )
    