clear all; close all; clc;

%% Parameter
r = 0.1133;                 % (% p.a.) discount rate
muy = 0.049; sy = 0.0675;   % (% p.a.) drift rate & volatility fare revenues
mux = 0.0725; sx = 0.2;     % (% p.a.) drift rate & volatility market prices after rail
rho = 0.8;                  % () correlation coefficient between y & x
I1 = 24.849;                % (€) one time construction cost infrastructure
c = 0.9;                    % (€ p.a.) operating costs infrastructre
I2 = 14.5;                  % (€) construction costs of property
K = 4.53;                   % (€) purchase costs of property development rights
p = 0.51;                   % (%) profit sharing rule p in (0,1)
epsilon = 15;               % penalty

% Plot Panel
[xst, ybar, xbar, ybare, xbare, ~, others] = thresholds(r,muy,sy,mux,sx,I1,c,I2,K,p,epsilon);
fprintf('Thresholds: \nx* = %.2f \nybar = %.2f \nybare = %.2f \nxbar = %.2f \nxbare = %.2f \n', xst, ybar, ybare, xbar, xbare)
y = linspace(0, ybar*1.02, 100)'; 
x = linspace(0, max(xst,xbar)*1.02, 100)';

[Vy, Vye, ~, ~] = arrayfun(@(y) boundaryV(r,muy,sy,mux,sx,I1,c,I2,K,p,epsilon,y,0), y);
[~, ~, Vx, Vxe] = arrayfun(@(x) boundaryV(r,muy,sy,mux,sx,I1,c,I2,K,p,epsilon,0,x), x);
[Vyybar, ~, ~, ~] = arrayfun(@(y) boundaryV(r,muy,sy,mux,sx,I1,c,I2,K,p,epsilon,y,0), ybar);
[~, Vyeybare, ~, ~] = arrayfun(@(y) boundaryV(r,muy,sy,mux,sx,I1,c,I2,K,p,epsilon,y,0), ybare);
[~, ~, Vxxbar, ~] = arrayfun(@(x) boundaryV(r,muy,sy,mux,sx,I1,c,I2,K,p,epsilon,0,x), xbar);
[~, ~, ~, Vxexbare] = arrayfun(@(x) boundaryV(r,muy,sy,mux,sx,I1,c,I2,K,p,epsilon,0,x), xbare);
phiy = arrayfun(@(y) obstacle(r,muy,sy,mux,sx,I1,c,I2,K,p,y,0), y);
phix = arrayfun(@(x) obstacle(r,muy,sy,mux,sx,I1,c,I2,K,p,0,x), x);
phixst = arrayfun(@(x) obstacle(r,muy,sy,mux,sx,I1,c,I2,K,p,0,x), xst);

figure(1)
subplot(1,2,1)
plot(y,Vy,'-k','Linewidth',1); hold on; plot(y,Vye,'--k','Linewidth',1); hold on; plot(y,phiy,'-','Color',[0.4 0.4 0.4],'Linewidth',1); hold on;
y1 = get(gca,'ylim'); x1 = get(gca,'xlim'); hold on;
plot([ybar ybar], [0 Vyybar], '-k'); hold on; plot([ybare ybare], [0 Vyeybare], '--k'); hold on;
plot(x1,[0,0],'-k','Linewidth',1); plot(x1,[-others.b -others.b],'-k');
set(gca,'FontSize',12,'FontUnits','points','FontName','times',...
'XColor','none','YColor',[0 0 0],'Linewidth',1,...
'box','off','XTick',[],'YTick',[0]); hold on;
axis tight;

subplot(1,2,2)
plot(x,Vx,'-k','Linewidth',1); hold on; plot(x,Vxe,'--k','Linewidth',1); hold on; plot(x,phix,'-','Color',[0.4 0.4 0.4],'Linewidth',1); hold on;
y1 = get(gca,'ylim'); x1 = get(gca,'xlim'); hold on;
xlabel('x');
plot([xbar xbar], [0 Vxxbar], '-k'); hold on; plot([xbare xbare], [0 Vxexbare], '--k'); hold on; plot([xst xst], [0 phixst], '-','Color',[0.4 0.4 0.4]); hold on;
plot(x1,[0,0],'-k','Linewidth',1); hold on; plot(x1,[-others.b -others.b],'-k'); xlim([0 x(end)]);
set(gca,'FontSize',12,'FontUnits','points','FontName','times',...
'XColor','none','YColor',[0 0 0],'Linewidth',1,'box','off','XTick',[],'YTick',[0]);
axis tight;

% Manual labeling is required for optimum positioning of the labels.

%% Define functions
function [payoff] = obstacle(r,muy,sy,mux,sx,I1,c,I2,K,p,y,x)
    
    % define a, alpha, b and beta
    a = exp(-(r-muy))/(r-muy);
    b = c/r*exp(-r)+I1;
    alpha = p*exp(-(r-mux));
    beta = p*I2-K;

    % case y = 0
    gamma = 1/2-mux/sx^2+sqrt((1/2-mux/sx^2)^2+2*r/sx^2);
    xst = gamma/(gamma-1)*(I2+K/(1-p))*exp((r-mux));
    
    % leader payoff
    payoff = a.*y-b+((alpha.*xst-beta).*(x./xst).^gamma).*(x<xst)+(alpha.*x-beta).*(x>=xst);
    
end


function [xst, ybar, xbar, ybare, xbare, roots, others] = thresholds(r,muy,sy,mux,sx,I1,c,I2,K,p,epsilon)
    
    % define a, alpha, b and beta
    a = exp(-(r-muy))/(r-muy);
    b = c/r*exp(-r)+I1;
    alpha = p*exp(-(r-mux));
    beta = p*I2-K;
    
    % case x = 0
    delta = 1/2-muy/sy^2+sqrt((1/2-muy/sy^2)^2+2*r/sy^2);
    ybar = delta/(delta-1)*b/a;

    % case y = 0
    gamma = 1/2-mux/sx^2+sqrt((1/2-mux/sx^2)^2+2*r/sx^2);
    xst = gamma/(gamma-1)*(I2+K/(1-p))*exp((r-mux));
    xbar = gamma/(gamma-1)*(b+beta)/alpha;
    
    % Boundary Cases dependent on epsilon
    % x = 0
    deltaep = 1/2-muy/sy^2+sqrt((1/2-muy/sy^2)^2+2*(r+1/epsilon)/sy^2);
    deltaen = 1/2-muy/sy^2-sqrt((1/2-muy/sy^2)^2+2*(r+1/epsilon)/sy^2);
    ybare = ( delta-deltaen-1/deltaep*2/(epsilon*sy^2) )/( delta-deltaen-1/(deltaep-1)*2/(epsilon*sy^2) )*b/a;
    % y = 0
    gammaep = 1/2-mux/sx^2+sqrt((1/2-mux/sx^2)^2+2*(r+1/epsilon)/sx^2);
    gammaen = 1/2-mux/sx^2-sqrt((1/2-mux/sx^2)^2+2*(r+1/epsilon)/sx^2);
    xbare = ( gamma-gammaen-1/gammaep*2/(epsilon*sx^2) )/( gamma-gammaen-1/(gammaep-1)*2/(epsilon*sx^2) )*(b+beta)/alpha;
    
    roots.delta = delta; roots.deltaep = deltaep; roots.deltaen = deltaen;
    roots.gamma = gamma; roots.gammaep = gammaep; roots.gammaen = gammaen;
    
    others.a = a; others.b = b; others.alpha = alpha; others.beta = beta;
    
end

function [Vy, Vye, Vx, Vxe] = boundaryV(r,muy,sy,mux,sx,I1,c,I2,K,p,epsilon,y,x)
    [xst, ybar, xbar, ybare, xbare, roots, others] = thresholds(r,muy,sy,mux,sx,I1,c,I2,K,p,epsilon);
    delta = roots.delta; deltaep = roots.deltaep; deltaen = roots.deltaen;
    gamma = roots.gamma; gammaep = roots.gammaep; gammaen = roots.gammaen;
    a = others.a; b = others.b; alpha = others.alpha; beta = others.beta;

    % case x = 0
    Vye =   ( (y/ybare)^delta*obstacle(r,muy,sy,mux,sx,I1,c,I2,K,p,ybare,0) ).*(y<ybare)+...
            (  1/(deltaep-deltaen)*2/(epsilon*sy^2)*y^deltaen*integral(@(xi) obstacle(r,muy,sy,mux,sx,I1,c,I2,K,p,xi,0)./xi.^(deltaen+1),ybare,y)+...
               1/(deltaep-deltaen)*2/(epsilon*sy^2)*y^deltaep*integral(@(xi) obstacle(r,muy,sy,mux,sx,I1,c,I2,K,p,xi,0)./xi.^(deltaep+1),y,inf)+...
               (y/ybare)^deltaen*( obstacle(r,muy,sy,mux,sx,I1,c,I2,K,p,ybare,0)-...
               1/(deltaep-deltaen)*2/(epsilon*sy^2)*ybare^deltaep*integral(@(xi) obstacle(r,muy,sy,mux,sx,I1,c,I2,K,p,xi,0)./xi.^(deltaep+1),ybare,inf) ) ).*(y>=ybare);
    Vy = (b/(delta-1)*(y/ybar)^delta).*(0<y&&y<ybar)+(a*y-b).*(y>=ybar);
           
    % case y = 0
    Vxe =   ( (x/xbare)^gamma*obstacle(r,muy,sy,mux,sx,I1,c,I2,K,p,0,xbare) ).*(x<xbare)+...
            (  1/(gammaep-gammaen)*2/(epsilon*sx^2)*x^gammaep*integral(@(xi) obstacle(r,muy,sy,mux,sx,I1,c,I2,K,p,0,xi)./xi.^(gammaep+1),x,inf)+...
               (x/xbare)^gammaen*obstacle(r,muy,sy,mux,sx,I1,c,I2,K,p,0,xbare)+...
               1/(gammaep-gammaen)*2/(epsilon*sx^2)*x^gammaen*(integral(@(xi) obstacle(r,muy,sy,mux,sx,I1,c,I2,K,p,0,xi)./xi.^(gammaen+1),xbare,x)-...
               xbare^(gammaep-gammaen)*integral(@(xi) obstacle(r,muy,sy,mux,sx,I1,c,I2,K,p,0,xi)./xi.^(gammaep+1),xbare,inf) ) ).*(x>=xbare); 
           
    Vx = ((alpha*max(xst,xbar)-(beta+b))*(x/max(xst,xbar))^gamma).*(0<x&&x<max(xst,xbar))...
        +(alpha*x-(beta+b)).*(x>=max(xst,xbar));

end