clc
clear all


%#################################
%Globalparameter initial festlegen
%#################################

% Symbolische Laufvariable für die Summen festlegen
syms t11 t12 t21 t22 k tau

t = [t11;t12;t21;t22;k;tau];

%Spätester Markteintrittszeitpunkt für die Company und den Competitor
%festlegen. Ein späterer Markteintrittszeitpunkt ist nicht möglich.
N = 3;

%Kalkulationszinssatz festlegen
r = 0.03;

%Maximale Anzahl an Szenarien festlegen
e_max = 2;

%Ergebnisvektor vorbelegen
Ergebnisse = zeros(e_max,40);

%Vektor der Innovations- und Immitationsparameter initial festlegen
p = zeros (e_max,2);
q = zeros (e_max,2);


%#############################
%globale Competitor Parameter
%#############################

%Innovator Profile Competitor
v_M = 1;

%Spillover Effekte Competitor
S_M = 0.25;

%##########################
%globale Company Parameter
%##########################

%Increasing/Decreasing Factor of Initial Sales Volume
a = 1.15;

%Decreasing Factor of FM-Probability (Wettbewerbskennwert für einen Markt)
b = 1.8;

%Prozentualer Anteil der Ausgaben für Marketing an den Gesamtausgaben
%M=0,05
%M=0,1
%M=0,2
M = 0.1;

% Marketing-Technologie Matrix definieren
A = [1.1 1.05 1.02;
    1.15 1.08 1.03;
    1.2 1.1 1.04];

%#########################################################
%Lokale Parameter für alle Szenarien generieren
%#########################################################

%Company Innovator profile festlegen und speichern
v_C = unifrnd(0,2,e_max,1);
Ergebnisse(:,8) = v_C;

%Knowledge Potentials der Company pro Innovation festlegen und speichern
K_m = unifrnd(0,6,e_max,2);
Ergebnisse(:,9) = K_m(:,1);
Ergebnisse(:,10) = K_m(:,2);

%Initial kumulierte Umsätze pro Markt während der Marktphase festlegen und speichern
m_0 = unifrnd(1e6,1e8,e_max,2);
Ergebnisse(:,11) = m_0(:,1);
Ergebnisse(:,12) = m_0(:,2);

%Anfangsauszahlungen der Company pro Innovation festlegen und speichern
A_0 = unifrnd(1e6,1e7,e_max,2);
Ergebnisse(:,13) = A_0(:,1);
Ergebnisse(:,14) = A_0(:,2);

% Initial Coefficient of Innovation pro Markt festlegen
p_0 = unifrnd(0.001,0.1,e_max,2);
Ergebnisse(:,15) = p_0(:,1);
Ergebnisse(:,16) = p_0(:,2);

% Initial Coefficient of Imitation pro Markt festlegen
q_0 = unifrnd(0.2,1.5,e_max,2);
Ergebnisse(:,17) = q_0(:,1);
Ergebnisse(:,18) = q_0(:,2);

%Einen Preis der Innovation initial auf 1 setzen um für den Umsatz die
%korrekte Einheit zu haben. Für das weitere Vorgehen nicht relevant.
P = 1;

%Spillover Effekte der Company für beide Innovationen und alle Szenarien festlegen
%Vektor vorbelegen
S_C = zeros(e_max,2);

for i = 1:e_max
    %für Innovation 1
    if (K_m(i,1) > 5 && K_m(i,1) <= 6)
        S_C(i,1) = 0.1;
    elseif (K_m(i,1) > 2 && K_m(i,1) <= 5)
        S_C(i,1) = 0.25;
    else
        S_C(i,1) = 0.6;
    end
    
    %für Innovation 2
    if (K_m(i,2) > 5 && K_m(i,2) <= 6)
        S_C(i,2) = 0.1;
    elseif (K_m(i,2) > 2 && K_m(i,2) <= 5)
        S_C(i,2) = 0.25;
    else
        S_C(i,2) = 0.6;
    end
end

%Spillover Effekte für beide Innovationen im Ergebnisvektor speichern
Ergebnisse(:,19) = S_C(:,1);
Ergebnisse(:,20) = S_C(:,2);


%#############################
%Schleife über alle Szenarien
%#############################

for n = 1 : e_max
   
    %Coefficient of Innovation und Imitation pro Innovation und Markt in Abhängigkeit der Marketingmaßnahmen festlegen
    %Innovation 1 in Markt 1
    if (M == 0.05 && K_m(n,1) > 0 && K_m(n,1) <= 2)
        p(n,1) = p_0(n,1) * A(1,1);
        q(n,1) = q_0(n,1) * A(1,1);
    elseif (M == 0.1 && K_m(n,1) > 0 && K_m(n,1) <= 2)
        p(n,1) = p_0(n,1) * A(2,1);
        q(n,1) = q_0(n,1) * A(2,1);
    elseif (M == 0.2 && K_m(n,1) > 0 && K_m(n,1) <= 2)
        p(n,1) = p_0(n,1) * A(3,1);
        q(n,1) = q_0(n,1) * A(3,1);
        
    elseif (M == 0.05 && K_m(n,1) > 2 && K_m(n,1) < 5)
        p(n,1) = p_0(n,1) * A(1,2);
        q(n,1) = q_0(n,1) * A(1,2);
    elseif (M == 0.1 && K_m(n,1) > 2 && K_m(n,1) < 5)
        p(n,1) = p_0(n,1) * A(2,2);
        q(n,1) = q_0(n,1) * A(2,2);
    elseif (M == 0.2 && K_m(n,1) > 2 && K_m(n,1) < 5)
        p(n,1) = p_0(n,1) * A(3,2);
        q(n,1) = q_0(n,1) * A(3,2);
        
    elseif (M == 0.05 && K_m(n,1) > 5 && K_m(n,1) <= 6)
        p(n,1) = p_0(n,1) * A(1,3);
        q(n,1) = q_0(n,1) * A(1,3);
    elseif (M == 0.1 && K_m(n,1) > 5 && K_m(n,1) <= 6)
        p(n,1) = p_0(n,1) * A(2,3);
        q(n,1) = q_0(n,1) * A(2,3);
    elseif (M == 0.2 && K_m(n,1) > 5 && K_m(n,1) <= 6)
        p(n,1) = p_0(n,1) * A(3,3);
        q(n,1) = q_0(n,1) * A(3,3);
    end
    
        %Innovation 2 in Markt 1
    if (M == 0.05 && K_m(n,2) > 0 && K_m(n,2) <= 2)
        p(n,1) = p_0(n,1) * A(1,1);
        q(n,1) = q_0(n,1) * A(1,1);
    elseif (M == 0.1 && K_m(n,2) > 0 && K_m(n,2) <= 2)
        p(n,1) = p_0(n,1) * A(2,1);
        q(n,1) = q_0(n,1) * A(2,1);
    elseif (M == 0.2 && K_m(n,2) > 0 && K_m(n,2) <= 2)
        p(n,1) = p_0(n,1) * A(3,1);
        q(n,1) = q_0(n,1) * A(3,1);
        
    elseif (M == 0.05 && K_m(n,2) > 2 && K_m(n,2) < 5)
        p(n,1) = p_0(n,1) * A(1,2);
        q(n,1) = q_0(n,1) * A(1,2);
    elseif (M == 0.1 && K_m(n,2) > 2 && K_m(n,2) < 5)
        p(n,1) = p_0(n,1) * A(2,2);
        q(n,1) = q_0(n,1) * A(2,2);
    elseif (M == 0.2 && K_m(n,2) > 2 && K_m(n,2) < 5)
        p(n,1) = p_0(n,1) * A(3,2);
        q(n,1) = q_0(n,1) * A(3,2);
        
    elseif (M == 0.05 && K_m(n,2) > 5 && K_m(n,2) <= 6)
        p(n,1) = p_0(n,1) * A(1,3);
        q(n,1) = q_0(n,1) * A(1,3);
    elseif (M == 0.1 && K_m(n,2) > 5 && K_m(n,2) <= 6)
        p(n,1) = p_0(n,1) * A(2,3);
        q(n,1) = q_0(n,1) * A(2,3);
    elseif (M == 0.2 && K_m(n,2) > 5 && K_m(n,2) <= 6)
        p(n,1) = p_0(n,1) * A(3,3);
        q(n,1) = q_0(n,1) * A(3,3);
    end
        
        %Innovation 1 in Markt 2
    if (M == 0.05 && K_m(n,1) > 0 && K_m(n,1) <= 2)
        p(n,2) = p_0(n,2) * A(1,1);
        q(n,2) = q_0(n,2) * A(1,1);
    elseif (M == 0.1 && K_m(n,1) > 0 && K_m(n,1) <= 2)
        p(n,2) = p_0(n,2) * A(2,1);
        q(n,2) = q_0(n,2) * A(2,1);
    elseif (M == 0.2 && K_m(n,1) > 0 && K_m(n,1) <= 2)
        p(n,2) = p_0(n,2) * A(3,1);
        q(n,2) = q_0(n,2) * A(3,1);
        
    elseif (M == 0.05 && K_m(n,1) > 2 && K_m(n,1) < 5)
        p(n,2) = p_0(n,2) * A(1,2);
        q(n,2) = q_0(n,2) * A(1,2);
    elseif (M == 0.1 && K_m(n,1) > 2 && K_m(n,1) < 5)
        p(n,2) = p_0(n,2) * A(2,2);
        q(n,2) = q_0(n,2) * A(2,2);
    elseif (M == 0.2 && K_m(n,1) > 2 && K_m(n,1) < 5)
        p(n,2) = p_0(n,2) * A(3,2);
        q(n,2) = q_0(n,2) * A(3,2);
        
    elseif (M == 0.05 && K_m(n,1) > 5 && K_m(n,1) <= 6)
        p(n,2) = p_0(n,2) * A(1,3);
        q(n,2) = q_0(n,2) * A(1,3);
    elseif (M == 0.1 && K_m(n,1) > 5 && K_m(n,1) <= 6)
        p(n,2) = p_0(n,2) * A(2,3);
        q(n,2) = q_0(n,2) * A(2,3);
    elseif (M == 0.2 && K_m(n,1) > 5 && K_m(n,1) <= 6)
        p(n,2) = p_0(n,2) * A(3,3);
        q(n,2) = q_0(n,2) * A(3,3);
    end
        
        %Innovation 2 in Markt 2
    if (M == 0.05 && K_m(n,2) > 0 && K_m(n,2) <= 2)
        p(n,2) = p_0(n,2) * A(1,1);
        q(n,2) = q_0(n,2) * A(1,1);
    elseif (M == 0.1 && K_m(n,2) > 0 && K_m(n,2) <= 2)
        p(n,2) = p_0(n,2) * A(2,1);
        q(n,2) = q_0(n,2) * A(2,1);
    elseif (M == 0.2 && K_m(n,2) > 0 && K_m(n,2) <= 2)
        p(n,2) = p_0(n,2) * A(3,1);
        q(n,2) = q_0(n,2) * A(3,1);
        
    elseif (M == 0.05 && K_m(n,2) > 2 && K_m(n,2) < 5)
        p(n,2) = p_0(n,2) * A(1,2);
        q(n,2) = q_0(n,2) * A(1,2);
    elseif (M == 0.1 && K_m(n,2) > 2 && K_m(n,2) < 5)
        p(n,2) = p_0(n,2) * A(2,2);
        q(n,2) = q_0(n,2) * A(2,2);
    elseif (M == 0.2 && K_m(n,2) > 2 && K_m(n,2) < 5)
        p(n,2) = p_0(n,2) * A(3,2);
        q(n,2) = q_0(n,2) * A(3,2);
        
    elseif (M == 0.05 && K_m(n,2) > 5 && K_m(n,2) <= 6)
        p(n,2) = p_0(n,2) * A(1,3);
        q(n,2) = q_0(n,2) * A(1,3);
    elseif (M == 0.1 && K_m(n,2) > 5 && K_m(n,2) <= 6)
        p(n,2) = p_0(n,2) * A(2,3);
        q(n,2) = q_0(n,2) * A(2,3);
    elseif (M == 0.2 && K_m(n,2) > 5 && K_m(n,2) <= 6)
        p(n,2) = p_0(n,2) * A(3,3);
        q(n,2) = q_0(n,2) * A(3,3);
    end
        
    
    %######################
    %Funktionen definieren
    %######################
    
    % Company Auszahlungen für R&D
    CF_R_11 = A_0(n,1) / 10 * K_m(n,1);
    CF_R_12 = A_0(n,1) / 10 * K_m(n,1);
    CF_R_21 = A_0(n,2) / 10 * K_m(n,2);
    CF_R_22 = A_0(n,2) / 10 * K_m(n,2);
    
    % Company Auszahlungen für Marketing zum Eintrittszeitpunkt
    CF_M_11 = M * CF_R_11 * t11;
    CF_M_12 = M * CF_R_12 * t12;
    CF_M_21 = M * CF_R_21 * t21;
    CF_M_22 = M * CF_R_22 * t22;
    
    % Verschiebungsparameter festlegen
    if (v_C(n,1) > 0 && v_C(n,1) < 1)
        z_C_11 = (1 - 0.5 * v_C(n,1)) * t11;
        z_C_12 = (1 - 0.5 * v_C(n,1)) * t12;
        z_C_21 = (1 - 0.5 * v_C(n,1)) * t21;
        z_C_22 = (1 - 0.5 * v_C(n,1)) * t22;
    else
        z_C_11 = 1 / (2 * v_C(n,1)) * t11;
        z_C_12 = 1 / (2 * v_C(n,1)) * t12;
        z_C_21 = 1 / (2 * v_C(n,1)) * t21;
        z_C_22 = 1 / (2 * v_C(n,1)) * t22;
    end
    
    % Company Wissen definieren
    K_C_11 = S_C(n,1) * K_m(n,1) - K_m(n,1) * (1 - S_C(n,1)) / (1 + exp(v_C(n,1) * z_C_11))...
        + K_m(n,1) * (1 - S_C(n,1)) / (1 + exp(-v_C(n,1) * (t11 - z_C_11)));
    
    K_C_12 = S_C(n,1) * K_m(n,1) - K_m(n,1) * (1 - S_C(n,1)) / (1 + exp(v_C(n,1) * z_C_12))...
        + K_m(n,1) * (1 - S_C(n,1)) / (1 + exp(-v_C(n,1) * (t12 - z_C_12)));
    
    K_C_21 = S_C(n,2) * K_m(n,2) - K_m(n,2) * (1 - S_C(n,2)) / (1 + exp(v_C(n,1) * z_C_21))...
        + K_m(n,2) * (1 - S_C(n,2)) / (1 + exp(-v_C(n,1) * (t21 - z_C_21)));
    
    K_C_22 = S_C(n,2) * K_m(n,2) - K_m(n,2) * (1 - S_C(n,2)) / (1 + exp(v_C(n,1) * z_C_22))...
        + K_m(n,2) * (1 - S_C(n,2)) / (1 + exp(-v_C(n,1) * (t22 - z_C_22)));
    

    % Competitor Wissen definieren
    K_M_11 = 0.25 * K_m(n,1) - 0.75 * K_m(n,1) / (1 + exp(0.5 * tau))...
        + 0.75 * K_m(n,1) / (1 + exp(-0.5 * tau));
    
    K_M_12 = 0.25 * K_m(n,1) - 0.75 * K_m(n,1) / (1 + exp(0.5 * tau))...
        + 0.75 * K_m(n,1) / (1 + exp(-0.5 * tau));
    
    K_M_21 = 0.25 * K_m(n,2) - 0.75 * K_m(n,2) / (1 + exp(0.5 * tau))...
        + 0.75 * K_m(n,2) / (1 + exp(-0.5 * tau));
    
    K_M_22 = 0.25 * K_m(n,2) - 0.75 * K_m(n,2) / (1 + exp(0.5 * tau))...
        + 0.75 * K_m(n,2) / (1 + exp(-0.5 * tau));
    
    % Company Sales Potential als First Mover definieren
    m_FM_C_11 = m_0(n,1) * a^K_C_11;
    m_FM_C_12 = m_0(n,2) * a^K_C_12;
    m_FM_C_21 = m_0(n,1) * a^K_C_21;
    m_FM_C_22 = m_0(n,2) * a^K_C_22;
    
    % Competitor Sales Potential als First Mover definieren
    m_FM_M_11 = m_0(n,1) * a^K_M_11;
    m_FM_M_12 = m_0(n,2) * a^K_M_12;
    m_FM_M_21 = m_0(n,1) * a^K_M_21;
    m_FM_M_22 = m_0(n,2) * a^K_M_22;
    
    % Company Sales Potential als Late Mover definieren
    m_LM_C_11 = m_0(n,1) * a^(K_C_11 - t11 + tau);
    m_LM_C_12 = m_0(n,2) * a^(K_C_12 - t12 + tau);
    m_LM_C_21 = m_0(n,1) * a^(K_C_21 - t21 + tau);
    m_LM_C_22 = m_0(n,2) * a^(K_C_22 - t22 + tau);
    
    
    % Erfolgswahrscheinlichkeiten für die Innovationen in den Märkten definieren
    Pr_S_11 = 1.1^(p(n,1) + q(n,1) + K_C_11) * t11 / ...
        (1.1^(p(n,1) + q(n,1) + K_m(n,1)) * t11 + 1.4^(-t11));
    
    Pr_S_12 = 1.1^(p(n,2) + q(n,2) + K_C_12) * t12 / ...
        (1.1^(p(n,2) + q(n,2) + K_m(n,1)) * t12 + 1.4^(-t12));
    
    Pr_S_21 = 1.1^(p(n,1) + q(n,1) + K_C_21) * t21 / ...
        (1.1^(p(n,1) + q(n,1) + K_m(n,2)) * t21 + 1.4^(-t21));
    
    Pr_S_22 = 1.1^(p(n,2) + q(n,2) + K_C_22) * t22 / ...
        (1.1^(p(n,2) + q(n,2) + K_m(n,2)) * t22 + 1.4^(-t22));
    

    % Scheiterwahrscheinlichkeit für die Innovationen in den Märkten
    % definieren
    Pr_F_11 = 1 - Pr_S_11;
    Pr_F_12 = 1 - Pr_S_12;
    Pr_F_21 = 1 - Pr_S_21;
    Pr_F_22 = 1 - Pr_S_22;
    
    % S-Kurve für Cashflows der Company und des Competitor pro Markt und Innovation definieren
    
    % S-Kurven der Company definieren
    S_C_11 = (p(n,1) + q(n,1))^2 / p(n,1) * ...
        exp(-(p(n,1) + q(n,1)) * (k - t11)) / (1 + (q(n,1) / p(n,1)) * ...
        exp(-(p(n,1) + q(n,1)) * (k - t11)));
    
    S_C_12 = (p(n,2) + q(n,2))^2 / p(n,2) * ...
        exp(-(p(n,2) + q(n,2)) * (k - t12)) / (1 + (q(n,2) / p(n,2)) * ...
        exp(-(p(n,2) + q(n,2)) * (k - t12)));
    
    S_C_21 = (p(n,1) + q(n,1))^2 / p(n,1) * ...
        exp(-(p(n,1) + q(n,1)) * (k - t21)) / (1 + (q(n,1) / p(n,1)) * ...
        exp(-(p(n,1) + q(n,1)) * (k - t21)));
    
    S_C_22 = (p(n,2) + q(n,2))^2 / p(n,2) * ...
        exp(-(p(n,2) + q(n,2)) * (k - t22)) / (1 + (q(n,2) / p(n,2)) * ...
        exp(-(p(n,2) + q(n,2)) * (k - t22)));
    
    % S-Kurven des Competitor definieren
    S_M_11 = (p(n,1) + q(n,1))^2 / p(n,1) * ...
        exp(-(p(n,1) + q(n,1)) * (k - tau)) / (1 + (q(n,1) / p(n,1)) * ...
        exp(-(p(n,1) + q(n,1)) * (k - tau)));
    
    S_M_12 = (p(n,2) + q(n,2))^2 / p(n,2) * ...
        exp(-(p(n,2) + q(n,2)) * (k - tau)) / (1 + (q(n,2) / p(n,2)) * ...
        exp(-(p(n,2) + q(n,2)) * (k - tau)));
    
    S_M_21 = (p(n,1) + q(n,1))^2 / p(n,1) * ...
        exp(-(p(n,1) + q(n,1)) * (k - tau)) / (1 + (q(n,1) / p(n,1)) * ...
        exp(-(p(n,1) + q(n,1)) * (k - tau)));
    
    S_M_22 = (p(n,2) + q(n,2))^2 / p(n,2) * ...
        exp(-(p(n,2) + q(n,2)) * (k - tau)) / (1 + (q(n,2) / p(n,2)) * ...
        exp(-(p(n,2) + q(n,2)) * (k - tau)));
    
    % Cashflow Funktionen der Company pro Markt und Innovation definieren
    % First Mover Cashflows definieren
    % Funktionen abhängig von Laufindex k (siehe Paper)
    CF_FM_C_11_1 = m_FM_C_11 * S_C_11;
    CF_FM_C_11_2 = m_FM_C_11 * S_C_11 * (1 - S_M_11);
    
    CF_FM_C_21_1 = m_FM_C_21 * S_C_21;
    CF_FM_C_21_2 = m_FM_C_21 * S_C_21 * (1 - S_M_21);
    
    CF_FM_C_12_1 = m_FM_C_12 * S_C_12;
    CF_FM_C_12_2 = m_FM_C_12 * S_C_12 * (1 - S_M_12);
    
    CF_FM_C_22_1 = m_FM_C_22 * S_C_22;
    CF_FM_C_22_2 = m_FM_C_22 * S_C_22 * (1 - S_M_22);
    
    % Late Mover Cashflows definieren
    % abhängig von k und tau (siehe Paper)
    CF_LM_C_11 = S_C_11 * (m_LM_C_11 + m_FM_M_11 * S_M_11);
    
    CF_LM_C_21 = S_C_21 * (m_LM_C_21 + m_FM_M_21 * S_M_21);
    
    CF_LM_C_12 = S_C_12 * (m_LM_C_12 + m_FM_M_12 * S_M_12);
    
    CF_LM_C_22 = S_C_22 * (m_LM_C_22 + m_FM_M_22 * S_M_22);
   
        % ###### First-Mover-Wahrscheinlichkeiten ######
    p_FM_11 = b^(1 - t11);
    p_FM_12 = b^(1 - t12);
    p_FM_21 = b^(1 - t21);
    p_FM_22 = b^(1 - t22);
    
    % ###### Late-Mover-Wahrscheinlichkeiten ######
    p_LM_11 = 1 - p_FM_11;
    p_LM_12 = 1 - p_FM_12;
    p_LM_21 = 1 - p_FM_21;
    p_LM_22 = 1 - p_FM_22;

    % Durchschnittliche Erwartungswerte der Company-NPV´s definieren
    % First-Mover NPV´s
    NPV_FM_C_11 = symsum(Pr_S_11 * (-A_0(n,1) + symsum(-CF_R_11 / (1 + r)^k,k,1,t11-1) + ...
        symsum(CF_FM_C_11_1 / (1 + r)^k,k,t11,tau-1) + ...
        symsum(CF_FM_C_11_2 / (1 + r)^k,k,tau,N)) + ...
        Pr_F_11 * (-A_0(n,1) + symsum(-CF_R_11 / (1 + r)^k,k,1,t11-1) + ...
        symsum(CF_FM_C_11_1 / (1 + r)^k,k,t11,t11)),tau,t11+1,N) / (N - t11);
    
    NPV_FM_C_12 = symsum(Pr_S_12 * (-A_0(n,1) + symsum(-CF_R_12 / (1 + r)^k, k,1,t12-1) + ...
        symsum(CF_FM_C_12_1 / (1 + r)^k,k,t12,tau-1) + ...
        symsum(CF_FM_C_12_2 / (1 + r)^k,k,tau,N)) + ...
        Pr_F_12 * (-A_0(n,1) + symsum(-CF_R_12 / (1 + r)^k, k,1,t12-1) + ...
        symsum(CF_FM_C_12_1 / (1 + r)^k,k,t12,t12)),tau,t12+1,N) / (N - t12);
    
    NPV_FM_C_21 = symsum(Pr_S_21 * (-A_0(n,2) + symsum(-CF_R_21 / (1 + r)^k, k,1,t21-1) + ...
        symsum(CF_FM_C_21_1 / (1 + r)^k, k,t21,tau-1) + ...
        symsum(CF_FM_C_21_2 / (1 + r)^k,k,tau,N)) + ...
        Pr_F_21 * (-A_0(n,2) + symsum(-CF_R_21 / (1 + r)^k, k,1,t21-1) + ...
        symsum(CF_FM_C_21_1 / (1 + r)^k,k,t21,t21)),tau,t21+1,N) / (N - t21);
    
    NPV_FM_C_22 = symsum(Pr_S_22 * (-A_0(n,2) + symsum(-CF_R_22 / (1 + r)^k, k,1,t22-1) + ...
        symsum(CF_FM_C_22_1 / (1 + r)^k, k,t22,tau-1) + ...
        symsum(CF_FM_C_22_2 / (1 + r)^k, k,tau,N)) + ...
        Pr_F_22 * (-A_0(n,2) + symsum(-CF_R_22 / (1 + r)^k, k,1,t22-1) + ...
        symsum(CF_FM_C_22_1 / (1 + r)^k,k,t22,t22)),tau,t22+1,N) / (N - t22);

    % Late-Mover NPV´s der Company
    NPV_LM_C_11 = symsum(Pr_S_11 * (-A_0(n,1) + symsum(-CF_R_11 / (1 + r)^k,k,1,t11-1) + ...
        symsum(CF_LM_C_11 / (1 + r)^k,k,t11,N)) + ...
        Pr_F_11 * (-A_0(n,1) + symsum(-CF_R_11 / (1 + r)^k, k,1,t11-1) + ...
        symsum(CF_LM_C_11 / (1 + r)^k,k,t11,t11)),tau,1,t11-1) / (N - t11);
    
    NPV_LM_C_12 = symsum(Pr_S_12 * (-A_0(n,1) + symsum(-CF_R_12 / (1 + r)^k,k,1,t12-1) + ...
        symsum(CF_LM_C_12 / (1 + r)^k,k,t12,N)) + ...
        Pr_F_12 * (-A_0(n,1) + symsum(-CF_R_12 / (1 + r)^k, k,1,t12-1) + ...
        symsum(CF_LM_C_12 / (1 + r)^k,k,t12,t12)),tau,1,t12-1) / (N - t12);
    
    NPV_LM_C_21 = symsum(Pr_S_21 * (-A_0(n,2) + symsum(-CF_R_21 / (1 + r)^k,k,1,t21-1) + ...
        symsum(CF_LM_C_21 / (1 + r)^k,k,t21,N)) + ...
        Pr_F_21 * (-A_0(n,2) + symsum(-CF_R_21 / (1 + r)^k, k,1,t21-1) + ...
        symsum(CF_LM_C_21 / (1 + r)^k,k,t21,t21)),tau,1,t21-1) / (N - t21);
    
    NPV_LM_C_22 = symsum(Pr_S_22 * (-A_0(n,2) + symsum(-CF_R_22 / (1 + r)^k,k,1,t22-1) + ...
        symsum(CF_LM_C_22 / (1 + r)^k,k,t22,N)) + ...
        Pr_F_22 * (-A_0(n,2) + symsum(-CF_R_22 / (1 + r)^k, k,1,t22-1) + ...
        symsum(CF_LM_C_22 / (1 + r)^k,k,t22,t22)),tau,1,t22-1) / (N - t22);
    

    % ###### Zielfunktion für die Optimierung ######
    Zielfunktion = -(p_FM_11 *  NPV_FM_C_11 + p_LM_11 * NPV_LM_C_11 + ...
        p_FM_12 *  NPV_FM_C_12 + p_LM_12 * NPV_LM_C_12 + ...
        p_FM_21 *  NPV_FM_C_21 + p_LM_21 * NPV_LM_C_21 + ...
        p_FM_22 *  NPV_FM_C_22 + p_LM_22 * NPV_LM_C_22);

 grad_Ziel = jacobian(Zielfunktion,[t11;t12;t21;t22])
 hess_Ziel = jacobian (grad_Ziel,[t11;t12;t21;t22])
 
fh = matlabFunction(Zielfunktion,grad_Ziel,hess_Ziel,'vars',{t});

options = optimoptions('fmincon', ...
    'GradObj', 'on', ...
    'Hessian', 'on', ...
    'Algorithm','trust-region-reflective', ...
    'Display','final');
    
% ###### optimize with fmincon ######
% ###### [X,FVAL,EXITFLAG,OUTPUT,LAMBDA,GRAD,HESSIAN] ######

% = fmincon(FUN,X0,A,B,Aeq,Beq,LB,UB,NONLCON,OPTIONS)
[tfinal,fval,exitflag,output] = fmincon(fh,[1;1;1;1],[],[],[],[],[1;1;1;1],[N;N;N;N],[],options)

end