Bai, Bensoussan, Briest, Chevalier-Roignant

Required packages

ClearAll

Needs["NDSolve`FEM`"]

Parameters used for all graphs

We introduce a small penalty parameter used, e.g., for penalization purposes:

ϵ=(10^(-5));

This parameter drives the size of the bounded domain in which the PDEs are solved

daggertobarratio=1.15;

Baseline case

We first consider the baseline case with parameter values provided in the Overleaf project.

Solutions of the nonlinear penalized ODEs and PDE

The below code constructs the three following functions for the baseline case:

The solution of the benchmark penalized ODE in y, namely (ψ^0):=(V^0)-(Φ^0). This serves as a boundary condition for the two-dimensional PDE.

The solution of the benchmark penalized ODE in x, namely (Subscript[ψ, 0]):=(Subscript[V, 0])-(Subscript[Φ, 0]). This serves as a boundary condition for the two-dimensional PDE.

The solution of the penalized PDE,

Clear[ψ0,ψ00,ψ,Γ];
Module[
{
(*PRIMARY PARAMETERS*)
(*Process parameters*)
σx=.2,
σy=.0675,
ρ=0.8,
μx=0.0725,
μy=0.049,
r=0.1133,
(*Profit parameters*)
I1=24.849,
c=.9,
I2=14.5,
K=4.53,
α=.51
},
(*SECONDARY PARAMETERS*)
(*Parameters for obstacle*)
a=(1/(r-μy));
b=(c/r)+I1;
β=α*I2-K;
γ=((-2 μx+((σx)^2)+2 (Sqrt[2 r ((σx)^2)+(((μx-(((σx)^2)/2)))^2)]))/(2 ((σx)^2)));
γtilde=((-2 μx+((σx)^2)-2 (Sqrt[2 r ((σx)^2)+(((μx-(((σx)^2)/2)))^2)]))/(2 ((σx)^2)));
δ=((-2 μy+((σy)^2)+2 (Sqrt[2 r ((σy)^2)+(((μy-(((σy)^2)/2)))^2)]))/(2 ((σy)^2)));
xstar=(γ/(γ-1))*(I2+(K/(1-α)));
yFB=(δ/(δ-1))*(b/a);
xFB=(γ/(γ-1))*((b+β)/α);
(*Parameters for benchmark problems*)
δe=((-2 μy+((σy)^2)+2 (Sqrt[-2 (-r-(1/ϵ)) ((σy)^2)+(((μy-(((σy)^2)/2)))^2)]))/(2 ((σy)^2)));
γe=((-2 μx+((σx)^2)+2 (Sqrt[-2 (-r-(1/ϵ)) ((σx)^2)+(((μx-(((σx)^2)/2)))^2)]))/(2 ((σx)^2)));
yFBe=((δ+ϵ*δe*(r-μy*δ)+ϵ*r*(δ-1))/(δ+ϵ*δe*(r-μy*δ)+ϵ*r*(δ-1)-1))*(b/a);
xFBe=((γ+ϵ*γe*(r-μx*γ)+ϵ*r*(γ-1))/(γ+ϵ*γe*(r-μx*γ)+ϵ*r*(γ-1)-1))*((b+β)/α);
(*Domain of PDE*)
ydagger=yFB*daggertobarratio;
xdagger=xFB*daggertobarratio;
xunderdagger=0;
Omega=Rectangle[{0,xunderdagger},{ydagger,xdagger}];

(*INTERMEDIATE FUNCTIONS*)
(*Obstacles*)
Phi00[x_]:=-b + ({{\[Piecewise], ({{((((x/xstar)))^γ)*(α*xstar-β), 0<x<=xstar}, {α*x-β, x>=xstar }})}});
Phi[y_,x_]:=a*y-b + ({{\[Piecewise], ({{((((x/xstar)))^γ)*(α*xstar-β), 0<x<=xstar}, {α*x-β, x>=xstar }})}});
(*The threshold below which the NPV is negative is*)
zeroNPV:=(y/.Solve[Phi[y,#]==0,y][[1]])&;
(*Define the function V^0(y)*)
V0[y_]:=
({{\[Piecewise], ({{b/(δ-1)*((((y/yFB)))^δ), 0<y<yFB}, {a*y-b, y>=yFB}})}});
(*Define the function V_0(x)*)
V00[x_]:=
({{\[Piecewise], ({{(α*Max[xstar,xFB]-(β+b))*((((x/(Max[xstar,xFB]))))^γ), 0<=x<Max[xstar,xFB]}, {α*x-(β+b), x>=Max[xstar,xFB]}})}});

(*DIFFERENTIAL OPERATORS*)
(*Differential operator of the ODE in y*)
L0[v_,y_]:=(1/2)((σy)^2) (y^2) D[v,{y,2}]+μy y D[v,y]-r v;
(*Differential operator of the ODE in x*)
L00[v_,x_]:=(1/2)((σx)^2) (x^2) D[v,{x,2}]+μx x D[v,x]-r v;
(*Differential operator L*)
L[v_,y_,x_]:=(1/2)((σy)^2) (y^2) D[v,{y,2}]+(1/2) ((σx)^2)(x^2) D[v,{x,2}]+ρ σy σx y x D[v,y,x]+μy y D[v,y]+μx x D[v,x]-r v;

(*BOUNDARY CONDITIONS OF THE BENCHMARK PENALIZED ODEs*)
(*Boundary conditions for the penalized ODE in y*)
bc0={
DirichletCondition[ψ[y]==b,y==0],DirichletCondition[ψ[y]==0,y==ydagger]
};
(*Boundary conditions for the penalized ODE in x*)
bc00={
DirichletCondition[ψ[x]==b,x==xunderdagger],DirichletCondition[ψ[x]==0,x==xdagger]
};
(*Boundary conditions for the penalized PDE*)
bc={
DirichletCondition[ψ[y,x]==V0[y]-(a y-b),x==xunderdagger],
DirichletCondition[ψ[y,x]==V00[x]-Phi00[x],y==0],
DirichletCondition[ψ[y,x]==0,x==xdagger],
DirichletCondition[ψ[y,x]==0,y==ydagger]};

(*DIFFERENTIAL EQUATIONS*)
(*Nonlinear ODE for the benchmark penalized problem in y*)
nonlinearODE0=L0[ψ[y],y]+(1/ϵ)Max[- ψ[y];0]==-L0[a*y-b,y];
(*Nonlinear ODE for the benchmark penalized problem in x*)
nonlinearODE00=L00[ψ[x],x]+(1/ϵ)Max[- ψ[x];0]==-L00[Phi00[x],x];
(*Nonlinear PDE for the two-dimensional penalized problem*)
nonlinearPDE=L[ψ[y,x],y,x]+(1/ϵ)Max[ -ψ[y,x];0]==-L[Phi[y,x],y,x];

(*NONLINEAR ODE/PDE SOLVERs*)
ψ0=NDSolveValue[
{nonlinearODE0,bc0},ψ,{y,0,ydagger},
Method->"FiniteElement"];
ψ00=NDSolveValue[
{nonlinearODE00,bc00},ψ,{x,xunderdagger,xdagger},
Method->"FiniteElement"];
V00e[x_]:=ψ00[x]+Phi00[x];
ψ=NDSolveValue[
{nonlinearPDE,bc},ψ,{y,x}∈Omega,
Method->"FiniteElement"
];
Γ:=y/.FindRoot[ψ[y,#]==0,{y,yFB-(yFB/xFB)*#}]&;
];

AbsoluteTiming[%]

Export["ψbaseline.wxf",ψ];
Export["xdagger.wxf",xdagger];
Export["ydagger.wxf",ydagger];

Plots of the solutions of the nonlinear penalized ODEs and PDE

We now plot these solutions (against the corresponding obstacles) in the baseline case

Benchmark penalized ODE in y

Plot[{a*y-b,ψ0[y]+(a*y-b),1+V0[y]},{y,0,ydagger},
PlotLegends->Placed[{"y↦ay-b","y↦(Subscript[V, ϵ]^0)(y)","y↦(V^0)(y)"},{Left,Top}],
Frame->True,
FrameTicks->{{{{-b,"-b"}},None},{{{yFBe,"(Subscript[(OverBar[y]), ϵ])"}},{{yFBe,"(OverBar[y])"}}}},
GridLines->{{yFB,yFB},{0}},
GridLinesStyle->Directive[Gray,Dashed]]

Benchmark penalized ODE in x

Plot[{Phi00[x],V00e[x],1+V00[x]},{x,0,xdagger},
PlotLegends->Placed[{"x↦(Subscript[Φ, 0])(x)","x↦(Subscript[V, 0  ϵ])(x)","x↦(Subscript[V, 0])(x)"},{Left,Top}],
Frame->True,
FrameTicks->{{{{-b,"-b"}},None},{{{xFBe,"(Subscript[(OverBar[x]), ϵ])"}},{{xFB,"(OverBar[x])"}}}},
GridLines->{{xFBe},{0}},
GridLinesStyle->Directive[Gray,Dashed]]

Benchmark penalized PDE

(*We want to show convergence for x very small*)
With[{xsmall=0,xlarge=10,ylarge=ydagger},
boundaries=Plot[Γ[x],{x,xsmall,xlarge},
FrameLabel->{"Property price x","Fare revenues y"},(*Add frame labels*)
PlotLegends->Placed[{"x↦(Subscript[Γ, ϵ])(x)","x↦Φ(·,x)^(-1)(0)"},{Right,Center}],
PlotRange->{{xsmall,xlarge},{0,ylarge}},(*Ensures that all relevant parts of the region are shown*)
GridLines->{{xFB,xstar},{Γ[xsmall]}},
Frame->True,
FrameTicks->{{{{Γ[xsmall],"(Subscript[Γ, ϵ])(0)"}},{{yFB,"(Subscript[(OverBar[y]), ϵ])"}}},{All,None}}
]
]

Γ[xFB]

(*We also check the behavior at (OverBar[x])*)
Γ[xFB]

(*We now plot the free boundary curve*)With[{xsmall=0,xlarge=300},
boundaries=Plot[Γ[x],{x,xsmall,xlarge},
FrameLabel->{"Property price x","Fare revenues y"},(*Add frame labels*)
PlotLegends->Placed[{"x↦(Subscript[Γ, ϵ])(x)"},{Right,Center}],
PlotRange->{{xsmall,xlarge},{0,Automatic}},(*Ensures that all relevant parts of the region are shown*)
GridLines->{{xFB,xstar},{yFB,Γ[xsmall],Γ[3]}},
Frame->True,
FrameTicks->{{{{Γ[xsmall],"(Subscript[Γ, ϵ])(0)"},{Γ[3],"(Subscript[Γ, ϵ])(3)"},0},{{yFB,"(Subscript[(OverBar[y]), ϵ])"}}},{{{xstar,"(SubStar[x])"},{xFB,"(Subscript[(OverBar[x]), ϵ])"}},None}}
]
]

(*We now show the result for this particular case*)
Figure2=With[{xsmall=3,xlarge=300,ylarge=ydagger},
boundaries=ContourPlot[{ψ[y,x]==(1/(10^4)),Phi[y,x]==0},{x,xsmall,xlarge},{y,0,ylarge},
FrameLabel->{"Property price x","Fare revenues y"},(*Add frame labels*)
PlotLegends->Placed[{"x↦(Subscript[Γ, ϵ])(x)","x↦Φ(·,x)^(-1)(0)"},After],
PlotRange->{{xsmall,xlarge},{0,ylarge}},(*Ensures that all relevant parts of the region are shown*)
GridLines->{{xFBe,xstar},{Γ[xsmall]}},
FrameTicks->{{{{Γ[xsmall],"(Subscript[(OverBar[y]), ϵ])"}},None},{{{xFB,"(Subscript[(OverBar[x]), ϵ])"},{xstar,"(SuperStar[x])"}},None}},
Frame->True,
ContourStyle->{Directive[Black,Thickness[0.005]],(*Style for \(\psi=1/10^4\)*)Directive[Black,Dashed,Thickness[0.005]]}
];
region1a=RegionPlot[x>xFBe&&y<Γ[xsmall],{x,xsmall,xlarge},{y,0,ylarge},PlotStyle->Directive[GrayLevel[1.0,0.5]],BoundaryStyle->None];
region1b=RegionPlot[x>xFBe&&y>Γ[xsmall],{x,xsmall,xlarge},{y,0,ylarge},
PlotStyle->Directive[GrayLevel[0.95,0.5]],BoundaryStyle->None];region2a=RegionPlot[x<xstar&&y>Γ[xsmall],{x,0,xlarge},{y,0,ylarge},
PlotStyle->Directive[GrayLevel[0.90,0.5]],BoundaryStyle->None];
region2b=RegionPlot[xstar<x<xFBe&&y>Γ[xsmall],{x,0,xlarge},{y,0,ylarge},PlotStyle->Directive[GrayLevel[0.85,0.5]],BoundaryStyle->None];
region3=RegionPlot[ψ[y,x]>(1/(10^4)),{x,xsmall,xlarge},{y,0,ylarge},
PlotStyle->Directive[GrayLevel[0.0,0.5]],BoundaryStyle->None
];
region4a=RegionPlot[ψ[y,x]<=(1/(10^4))&&x<=xstar&&y<Γ[xsmall],{x,xsmall,xlarge},{y,0,ylarge},PlotStyle->Directive[GrayLevel[0.8,0.5]],BoundaryStyle->None];
region4b=RegionPlot[ψ[y,x]<=(1/(10^4))&&xstar<x<xFB&&y<Γ[xsmall],{x,xsmall,xlarge},{y,0,ylarge},PlotStyle->Directive[GrayLevel[0.8,0.5]],BoundaryStyle->None];
Show[
{boundaries,
region1a,region1b,region2a,region2b,region3,region4a,region4b,
Graphics[{{Text[Style["1a",Italic,Larger],{xFBe+25,((0+yFBe)/2)}],Text[Style["2a",Italic,Larger],{((0+xstar)/2),Γ[xsmall]+((ydagger-Γ[xsmall])/2)}],Text[Style["2b",Italic,Larger],{((xstar+xFBe)/2),Γ[xsmall]+((ydagger-Γ[xsmall])/2)}],Text[Style["1b/2c",Italic,Larger],{xFBe+25,Γ[xsmall]+((ydagger-Γ[xsmall])/2)}],Text[Style["3",Italic,Larger],{((0+xstar)/2),((0+yFBe)/2)}],Text[Style["4a",Italic,Larger],{((0+4*xstar)/5),Γ[xsmall]-.4}],Text[Style["4b",Italic,Larger],{((xstar+xFBe)/2),((0+yFBe)/2)}],
} }]
},
AspectRatio->2/3
]
];
Export["Figure2.png",Figure2];SystemOpen["Figure2.png"]

Comparative statics

aspectratio=(13/9);
xsmall=3;
xlarge=250;
ylarge=4;

Panel A -- Growth rate (Subscript[μ, x])

(Subscript[μ, x])=0.0000 (low)

Clear[ψ,xdagger,ydagger];
Module[
{
(*PRIMARY PARAMETERS*)
(*Process parameters*)
σx=.2,
σy=.0675,
ρ=0.8,
μx=0.0000, (*CHANGED PARAMETER*)
μy=0.049,
r=0.1133,
(*Profit parameters*)
I1=24.849,
c=.9,
I2=14.5,
K=4.53,
α=.51
},
(*SECONDARY PARAMETERS*)
(*Parameters for obstacle*)
a=(1/(r-μy));
b=(c/r)+I1;
β=α*I2-K;
γ=((-2 μx+((σx)^2)+2 (Sqrt[2 r ((σx)^2)+(((μx-(((σx)^2)/2)))^2)]))/(2 ((σx)^2)));
γtilde=((-2 μx+((σx)^2)-2 (Sqrt[2 r ((σx)^2)+(((μx-(((σx)^2)/2)))^2)]))/(2 ((σx)^2)));
δ=((-2 μy+((σy)^2)+2 (Sqrt[2 r ((σy)^2)+(((μy-(((σy)^2)/2)))^2)]))/(2 ((σy)^2)));
xstar=(γ/(γ-1))*(I2+(K/(1-α)));
yFB=(δ/(δ-1))*(b/a);
xFB=(γ/(γ-1))*((b+β)/α);
(*Parameters for benchmark problems*)
δe=((-2 μy+((σy)^2)+2 (Sqrt[-2 (-r-(1/ϵ)) ((σy)^2)+(((μy-(((σy)^2)/2)))^2)]))/(2 ((σy)^2)));
γe=((-2 μx+((σx)^2)+2 (Sqrt[-2 (-r-(1/ϵ)) ((σx)^2)+(((μx-(((σx)^2)/2)))^2)]))/(2 ((σx)^2)));
yFBe=((δ+ϵ*δe*(r-μy*δ)+ϵ*r*(δ-1))/(δ+ϵ*δe*(r-μy*δ)+ϵ*r*(δ-1)-1))*(b/a);
xFBe=((γ+ϵ*γe*(r-μx*γ)+ϵ*r*(γ-1))/(γ+ϵ*γe*(r-μx*γ)+ϵ*r*(γ-1)-1))*((b+β)/α);
(*Domain of PDE*)
ydagger=yFBe*daggertobarratio;
xdagger=xFBe*daggertobarratio;
Omega=Rectangle[{0,0},{ydagger,xdagger}];

(*INTERMEDIATE FUNCTIONS*)
(*Obstacles*)
Phi00[x_]:=-b + ({{\[Piecewise], ({{((((x/xstar)))^γ)*(α*xstar-β), 0<x<=xstar}, {α*x-β, x>=xstar }})}});
Phi[y_,x_]:=a*y-b + ({{\[Piecewise], ({{((((x/xstar)))^γ)*(α*xstar-β), 0<x<=xstar}, {α*x-β, x>=xstar }})}});
(*Define the function V^0(y)*)
V0[y_]:=
({{\[Piecewise], ({{b/(δ-1)*((((y/yFB)))^δ), 0<y<yFB}, {a*y-b, y>=yFB}})}});
(*Define the function V_0(x)*)
V00[x_]:=
({{\[Piecewise], ({{(α*Max[xstar,xFB]-(β+b))*((((x/(Max[xstar,xFB]))))^γ), 0<=x<Max[xstar,xFB]}, {α*x-(β+b), x>=Max[xstar,xFB]}})}});

(*DIFFERENTIAL OPERATORS*)
(*Differential operator of the ODE in y*)
L0[v_,y_]:=(1/2)((σy)^2) (y^2) D[v,{y,2}]+μy y D[v,y]-r v;
(*Differential operator of the ODE in x*)
L00[v_,x_]:=(1/2)((σx)^2) (x^2) D[v,{x,2}]+μx x D[v,x]-r v;
(*Differential operator L*)
L[v_,y_,x_]:=(1/2)((σy)^2) (y^2) D[v,{y,2}]+(1/2) ((σx)^2)(x^2) D[v,{x,2}]+ρ σy σx y x D[v,y,x]+μy y D[v,y]+μx x D[v,x]-r v;

(*BOUNDARY CONDITIONS OF THE BENCHMARK PENALIZED ODEs*)
(*Boundary conditions for the penalized ODE in y*)
bc0={
DirichletCondition[ψ[y]==b,y==0],DirichletCondition[ψ[y]==0,y==ydagger]
};
(*Boundary conditions for the penalized ODE in x*)
bc00={
DirichletCondition[ψ[x]==b,x==0],DirichletCondition[ψ[x]==0,x==xdagger]
};
(*Boundary conditions for the penalized PDE*)
bc={
DirichletCondition[ψ[y,x]==V0[y]-(a y-b),x==0],
DirichletCondition[ψ[y,x]==V00[x]-Phi00[x],y==0],
DirichletCondition[ψ[y,x]==0,x==xdagger],
DirichletCondition[ψ[y,x]==0,y==ydagger]};

(*DIFFERENTIAL EQUATIONS*)
(*Nonlinear ODE for the benchmark penalized problem in y*)
nonlinearODE0=L0[ψ[y],y]+(1/ϵ)Max[- ψ[y];0]==-L0[a*y-b,y];
(*Nonlinear ODE for the benchmark penalized problem in x*)
nonlinearODE00=L00[ψ[x],x]+(1/ϵ)Max[- ψ[x];0]==-L00[Phi00[x],x];
(*Nonlinear PDE for the two-dimensional penalized problem*)
nonlinearPDE=L[ψ[y,x],y,x]+(1/ϵ)Max[ -ψ[y,x];0]==-L[Phi[y,x],y,x];

(*NONLINEAR ODE/PDE SOLVERs*)
ψ0=NDSolveValue[
{nonlinearODE0,bc0},ψ,{y,0,ydagger},
Method->"FiniteElement"];
ψ00=NDSolveValue[
{nonlinearODE00,bc00},ψ,{x,0,xdagger},
Method->"FiniteElement"];
V00e[x_]:=ψ00[x]+Phi00[x];
ψ=NDSolveValue[
{nonlinearPDE,bc},ψ,{y,x}∈Omega,
Method->"FiniteElement"
];
];
Export["ψμx0.wxf",ψ];

(Subscript[μ, x])=0.0500 (medium)

Clear[ψ,xdagger,ydagger];
Module[
{
(*PRIMARY PARAMETERS*)
(*Process parameters*)
σx=.2,
σy=.0675,
ρ=0.8,
μx=0.0500, (*CHANGED PARAMETER*)
μy=0.049,
r=0.1133,
(*Profit parameters*)
I1=24.849,
c=.9,
I2=14.5,
K=4.53,
α=.51
},
(*SECONDARY PARAMETERS*)
(*Parameters for obstacle*)
a=(1/(r-μy));
b=(c/r)+I1;
β=α*I2-K;
γ=((-2 μx+((σx)^2)+2 (Sqrt[2 r ((σx)^2)+(((μx-(((σx)^2)/2)))^2)]))/(2 ((σx)^2)));
γtilde=((-2 μx+((σx)^2)-2 (Sqrt[2 r ((σx)^2)+(((μx-(((σx)^2)/2)))^2)]))/(2 ((σx)^2)));
δ=((-2 μy+((σy)^2)+2 (Sqrt[2 r ((σy)^2)+(((μy-(((σy)^2)/2)))^2)]))/(2 ((σy)^2)));
xstar=(γ/(γ-1))*(I2+(K/(1-α)));
yFB=(δ/(δ-1))*(b/a);
xFB=(γ/(γ-1))*((b+β)/α);
(*Parameters for benchmark problems*)
δe=((-2 μy+((σy)^2)+2 (Sqrt[-2 (-r-(1/ϵ)) ((σy)^2)+(((μy-(((σy)^2)/2)))^2)]))/(2 ((σy)^2)));
γe=((-2 μx+((σx)^2)+2 (Sqrt[-2 (-r-(1/ϵ)) ((σx)^2)+(((μx-(((σx)^2)/2)))^2)]))/(2 ((σx)^2)));
yFBe=((δ+ϵ*δe*(r-μy*δ)+ϵ*r*(δ-1))/(δ+ϵ*δe*(r-μy*δ)+ϵ*r*(δ-1)-1))*(b/a);
xFBe=((γ+ϵ*γe*(r-μx*γ)+ϵ*r*(γ-1))/(γ+ϵ*γe*(r-μx*γ)+ϵ*r*(γ-1)-1))*((b+β)/α);
(*Domain of PDE*)
ydagger=yFBe*daggertobarratio;
xdagger=xFBe*daggertobarratio;
Omega=Rectangle[{0,0},{ydagger,xdagger}];

(*INTERMEDIATE FUNCTIONS*)
(*Obstacles*)
Phi00[x_]:=-b + ({{\[Piecewise], ({{((((x/xstar)))^γ)*(α*xstar-β), 0<x<=xstar}, {α*x-β, x>=xstar }})}});
Phi[y_,x_]:=a*y-b + ({{\[Piecewise], ({{((((x/xstar)))^γ)*(α*xstar-β), 0<x<=xstar}, {α*x-β, x>=xstar }})}});
(*Define the function V^0(y)*)
V0[y_]:=
({{\[Piecewise], ({{b/(δ-1)*((((y/yFB)))^δ), 0<y<yFB}, {a*y-b, y>=yFB}})}});
(*Define the function V_0(x)*)
V00[x_]:=
({{\[Piecewise], ({{(α*Max[xstar,xFB]-(β+b))*((((x/(Max[xstar,xFB]))))^γ), 0<=x<Max[xstar,xFB]}, {α*x-(β+b), x>=Max[xstar,xFB]}})}});

(*DIFFERENTIAL OPERATORS*)
(*Differential operator of the ODE in y*)
L0[v_,y_]:=(1/2)((σy)^2) (y^2) D[v,{y,2}]+μy y D[v,y]-r v;
(*Differential operator of the ODE in x*)
L00[v_,x_]:=(1/2)((σx)^2) (x^2) D[v,{x,2}]+μx x D[v,x]-r v;
(*Differential operator L*)
L[v_,y_,x_]:=(1/2)((σy)^2) (y^2) D[v,{y,2}]+(1/2) ((σx)^2)(x^2) D[v,{x,2}]+ρ σy σx y x D[v,y,x]+μy y D[v,y]+μx x D[v,x]-r v;

(*BOUNDARY CONDITIONS OF THE BENCHMARK PENALIZED ODEs*)
(*Boundary conditions for the penalized ODE in y*)
bc0={
DirichletCondition[ψ[y]==b,y==0],DirichletCondition[ψ[y]==0,y==ydagger]
};
(*Boundary conditions for the penalized ODE in x*)
bc00={
DirichletCondition[ψ[x]==b,x==0],DirichletCondition[ψ[x]==0,x==xdagger]
};
(*Boundary conditions for the penalized PDE*)
bc={
DirichletCondition[ψ[y,x]==V0[y]-(a y-b),x==0],
DirichletCondition[ψ[y,x]==V00[x]-Phi00[x],y==0],
DirichletCondition[ψ[y,x]==0,x==xdagger],
DirichletCondition[ψ[y,x]==0,y==ydagger]};

(*DIFFERENTIAL EQUATIONS*)
(*Nonlinear ODE for the benchmark penalized problem in y*)
nonlinearODE0=L0[ψ[y],y]+(1/ϵ)Max[- ψ[y];0]==-L0[a*y-b,y];
(*Nonlinear ODE for the benchmark penalized problem in x*)
nonlinearODE00=L00[ψ[x],x]+(1/ϵ)Max[- ψ[x];0]==-L00[Phi00[x],x];
(*Nonlinear PDE for the two-dimensional penalized problem*)
nonlinearPDE=L[ψ[y,x],y,x]+(1/ϵ)Max[ -ψ[y,x];0]==-L[Phi[y,x],y,x];

(*NONLINEAR ODE/PDE SOLVERs*)
ψ0=NDSolveValue[
{nonlinearODE0,bc0},ψ,{y,0,ydagger},
Method->"FiniteElement"];
ψ00=NDSolveValue[
{nonlinearODE00,bc00},ψ,{x,0,xdagger},
Method->"FiniteElement"];
V00e[x_]:=ψ00[x]+Phi00[x];
ψ=NDSolveValue[
{nonlinearPDE,bc},ψ,{y,x}∈Omega,
Method->"FiniteElement"
];
];
Export["ψμx1.wxf",ψ];

Comparative statics graph for (Subscript[μ, x])

ψbaseline=Import["ψbaseline.wxf"];
xdagger=Import["xdagger.wxf"];
ydagger=Import["ydagger.wxf"];
ψ0=Import["ψμx0.wxf"];
ψ1=Import["ψμx1.wxf"];
Fig4a1=ContourPlot[ψ0[y,x]==(1/(10^4)),{x,xsmall,xlarge},{y,0,ydagger},
ContourStyle->{Directive[Black,Dotted,Thickness[0.005]]},
Frame->{True,True,False,False},
PlotRange->{{0,xlarge},{0,ylarge}},
FrameTicks->{{{1,2,3,4},None},{{0,50,100,150,200,250},None}},FrameLabel->{"Property price x","Fare revenues y","",""}
];
Fig4a2=ContourPlot[ψ1[y,x]==(1/(10^4)),{x,2,xdagger},{y,0,ydagger},Frame->True,
ContourStyle->{Directive[Black,Dashed,Thickness[0.005]]}];
Fig4a3=ContourPlot[ψbaseline[y,x]==(1/(10^4)),{x,xsmall,xdagger},{y,0,ydagger},ContourStyle->{Directive[Black,Plain,Thickness[0.005]]}];
Figure4a=Legended[
Show[{Fig4a1,Fig4a2,Fig4a3},AspectRatio->aspectratio],Placed[Framed@LineLegend[{Directive[Black,Dotted,Thickness[0.005]],Directive[Black,Dashed,Thickness[0.005]],Directive[Black,Plain,Thickness[0.005]]},{"(Subscript[μ, x])=.0000","(Subscript[μ, x])=.0500","(Subscript[μ, x])=0.0725"},LegendMarkerSize->15,Spacings->.2,LegendLabel->"Curve (Subscript[Γ, ϵ])(·) for"],{Center,Top}]];
Export["Figure4a.png",Figure4a];SystemOpen["Figure4a.png"]

Panel B -- Volatility (Subscript[σ, x])

(Subscript[σ, x])=0.15 (low)

Clear[ψ,xdagger,ydagger];
Module[
{
(*PRIMARY PARAMETERS*)
(*Process parameters*)
σx=.15, (*CHANGED PARAMETER*)
σy=.0675,
ρ=0.8,
μx=0.0725,
μy=0.049,
r=0.1133,
(*Profit parameters*)
I1=24.849,
c=.9,
I2=14.5,
K=4.53,
α=.51
},
(*SECONDARY PARAMETERS*)
(*Parameters for obstacle*)
a=(1/(r-μy));
b=(c/r)+I1;
β=α*I2-K;
γ=((-2 μx+((σx)^2)+2 (Sqrt[2 r ((σx)^2)+(((μx-(((σx)^2)/2)))^2)]))/(2 ((σx)^2)));
γtilde=((-2 μx+((σx)^2)-2 (Sqrt[2 r ((σx)^2)+(((μx-(((σx)^2)/2)))^2)]))/(2 ((σx)^2)));
δ=((-2 μy+((σy)^2)+2 (Sqrt[2 r ((σy)^2)+(((μy-(((σy)^2)/2)))^2)]))/(2 ((σy)^2)));
xstar=(γ/(γ-1))*(I2+(K/(1-α)));
yFB=(δ/(δ-1))*(b/a);
xFB=(γ/(γ-1))*((b+β)/α);
(*Parameters for benchmark problems*)
δe=((-2 μy+((σy)^2)+2 (Sqrt[-2 (-r-(1/ϵ)) ((σy)^2)+(((μy-(((σy)^2)/2)))^2)]))/(2 ((σy)^2)));
γe=((-2 μx+((σx)^2)+2 (Sqrt[-2 (-r-(1/ϵ)) ((σx)^2)+(((μx-(((σx)^2)/2)))^2)]))/(2 ((σx)^2)));
yFBe=((δ+ϵ*δe*(r-μy*δ)+ϵ*r*(δ-1))/(δ+ϵ*δe*(r-μy*δ)+ϵ*r*(δ-1)-1))*(b/a);
xFBe=((γ+ϵ*γe*(r-μx*γ)+ϵ*r*(γ-1))/(γ+ϵ*γe*(r-μx*γ)+ϵ*r*(γ-1)-1))*((b+β)/α);
(*Domain of PDE*)
ydagger=yFBe*daggertobarratio;
xdagger=xFBe*daggertobarratio;
Omega=Rectangle[{0,0},{ydagger,xdagger}];

(*INTERMEDIATE FUNCTIONS*)
(*Obstacles*)
Phi00[x_]:=-b + ({{\[Piecewise], ({{((((x/xstar)))^γ)*(α*xstar-β), 0<x<=xstar}, {α*x-β, x>=xstar }})}});
Phi[y_,x_]:=a*y-b + ({{\[Piecewise], ({{((((x/xstar)))^γ)*(α*xstar-β), 0<x<=xstar}, {α*x-β, x>=xstar }})}});
(*Define the function V^0(y)*)
V0[y_]:=
({{\[Piecewise], ({{b/(δ-1)*((((y/yFB)))^δ), 0<y<yFB}, {a*y-b, y>=yFB}})}});
(*Define the function V_0(x)*)
V00[x_]:=
({{\[Piecewise], ({{(α*Max[xstar,xFB]-(β+b))*((((x/(Max[xstar,xFB]))))^γ), 0<=x<Max[xstar,xFB]}, {α*x-(β+b), x>=Max[xstar,xFB]}})}});

(*DIFFERENTIAL OPERATORS*)
(*Differential operator of the ODE in y*)
L0[v_,y_]:=(1/2)((σy)^2) (y^2) D[v,{y,2}]+μy y D[v,y]-r v;
(*Differential operator of the ODE in x*)
L00[v_,x_]:=(1/2)((σx)^2) (x^2) D[v,{x,2}]+μx x D[v,x]-r v;
(*Differential operator L*)
L[v_,y_,x_]:=(1/2)((σy)^2) (y^2) D[v,{y,2}]+(1/2) ((σx)^2)(x^2) D[v,{x,2}]+ρ σy σx y x D[v,y,x]+μy y D[v,y]+μx x D[v,x]-r v;

(*BOUNDARY CONDITIONS OF THE BENCHMARK PENALIZED ODEs*)
(*Boundary conditions for the penalized ODE in y*)
bc0={
DirichletCondition[ψ[y]==b,y==0],DirichletCondition[ψ[y]==0,y==ydagger]
};
(*Boundary conditions for the penalized ODE in x*)
bc00={
DirichletCondition[ψ[x]==b,x==0],DirichletCondition[ψ[x]==0,x==xdagger]
};
(*Boundary conditions for the penalized PDE*)
bc={
DirichletCondition[ψ[y,x]==V0[y]-(a y-b),x==0],
DirichletCondition[ψ[y,x]==V00[x]-Phi00[x],y==0],
DirichletCondition[ψ[y,x]==0,x==xdagger],
DirichletCondition[ψ[y,x]==0,y==ydagger]};

(*DIFFERENTIAL EQUATIONS*)
(*Nonlinear ODE for the benchmark penalized problem in y*)
nonlinearODE0=L0[ψ[y],y]+(1/ϵ)Max[- ψ[y];0]==-L0[a*y-b,y];
(*Nonlinear ODE for the benchmark penalized problem in x*)
nonlinearODE00=L00[ψ[x],x]+(1/ϵ)Max[- ψ[x];0]==-L00[Phi00[x],x];
(*Nonlinear PDE for the two-dimensional penalized problem*)
nonlinearPDE=L[ψ[y,x],y,x]+(1/ϵ)Max[ -ψ[y,x];0]==-L[Phi[y,x],y,x];

(*NONLINEAR ODE/PDE SOLVERs*)
ψ0=NDSolveValue[
{nonlinearODE0,bc0},ψ,{y,0,ydagger},
Method->"FiniteElement"];
ψ00=NDSolveValue[
{nonlinearODE00,bc00},ψ,{x,0,xdagger},
Method->"FiniteElement"];
V00e[x_]:=ψ00[x]+Phi00[x];
ψ=NDSolveValue[
{nonlinearPDE,bc},ψ,{y,x}∈Omega,
Method->"FiniteElement"
];
];
Export["ψσx0.wxf",ψ];

(Subscript[σ, x])=0.25 (high)

Clear[ψ,xdagger,ydagger];
Module[
{
(*PRIMARY PARAMETERS*)
(*Process parameters*)
σx=.25, (*CHANGED PARAMETER*)
σy=.0675,
ρ=0.8,
μx=0.0725,
μy=0.049,
r=0.1133,
(*Profit parameters*)
I1=24.849,
c=.9,
I2=14.5,
K=4.53,
α=.51
},
(*SECONDARY PARAMETERS*)
(*Parameters for obstacle*)
a=(1/(r-μy));
b=(c/r)+I1;
β=α*I2-K;
γ=((-2 μx+((σx)^2)+2 (Sqrt[2 r ((σx)^2)+(((μx-(((σx)^2)/2)))^2)]))/(2 ((σx)^2)));
γtilde=((-2 μx+((σx)^2)-2 (Sqrt[2 r ((σx)^2)+(((μx-(((σx)^2)/2)))^2)]))/(2 ((σx)^2)));
δ=((-2 μy+((σy)^2)+2 (Sqrt[2 r ((σy)^2)+(((μy-(((σy)^2)/2)))^2)]))/(2 ((σy)^2)));
xstar=(γ/(γ-1))*(I2+(K/(1-α)));
yFB=(δ/(δ-1))*(b/a);
xFB=(γ/(γ-1))*((b+β)/α);
(*Parameters for benchmark problems*)
δe=((-2 μy+((σy)^2)+2 (Sqrt[-2 (-r-(1/ϵ)) ((σy)^2)+(((μy-(((σy)^2)/2)))^2)]))/(2 ((σy)^2)));
γe=((-2 μx+((σx)^2)+2 (Sqrt[-2 (-r-(1/ϵ)) ((σx)^2)+(((μx-(((σx)^2)/2)))^2)]))/(2 ((σx)^2)));
yFBe=((δ+ϵ*δe*(r-μy*δ)+ϵ*r*(δ-1))/(δ+ϵ*δe*(r-μy*δ)+ϵ*r*(δ-1)-1))*(b/a);
xFBe=((γ+ϵ*γe*(r-μx*γ)+ϵ*r*(γ-1))/(γ+ϵ*γe*(r-μx*γ)+ϵ*r*(γ-1)-1))*((b+β)/α);
(*Domain of PDE*)
ydagger=yFBe*daggertobarratio;
xdagger=xFBe*daggertobarratio;
Omega=Rectangle[{0,0},{ydagger,xdagger}];

(*INTERMEDIATE FUNCTIONS*)
(*Obstacles*)
Phi00[x_]:=-b + ({{\[Piecewise], ({{((((x/xstar)))^γ)*(α*xstar-β), 0<x<=xstar}, {α*x-β, x>=xstar }})}});
Phi[y_,x_]:=a*y-b + ({{\[Piecewise], ({{((((x/xstar)))^γ)*(α*xstar-β), 0<x<=xstar}, {α*x-β, x>=xstar }})}});
(*Define the function V^0(y)*)
V0[y_]:=
({{\[Piecewise], ({{b/(δ-1)*((((y/yFB)))^δ), 0<y<yFB}, {a*y-b, y>=yFB}})}});
(*Define the function V_0(x)*)
V00[x_]:=
({{\[Piecewise], ({{(α*Max[xstar,xFB]-(β+b))*((((x/(Max[xstar,xFB]))))^γ), 0<=x<Max[xstar,xFB]}, {α*x-(β+b), x>=Max[xstar,xFB]}})}});

(*DIFFERENTIAL OPERATORS*)
(*Differential operator of the ODE in y*)
L0[v_,y_]:=(1/2)((σy)^2) (y^2) D[v,{y,2}]+μy y D[v,y]-r v;
(*Differential operator of the ODE in x*)
L00[v_,x_]:=(1/2)((σx)^2) (x^2) D[v,{x,2}]+μx x D[v,x]-r v;
(*Differential operator L*)
L[v_,y_,x_]:=(1/2)((σy)^2) (y^2) D[v,{y,2}]+(1/2) ((σx)^2)(x^2) D[v,{x,2}]+ρ σy σx y x D[v,y,x]+μy y D[v,y]+μx x D[v,x]-r v;

(*BOUNDARY CONDITIONS OF THE BENCHMARK PENALIZED ODEs*)
(*Boundary conditions for the penalized ODE in y*)
bc0={
DirichletCondition[ψ[y]==b,y==0],DirichletCondition[ψ[y]==0,y==ydagger]
};
(*Boundary conditions for the penalized ODE in x*)
bc00={
DirichletCondition[ψ[x]==b,x==0],DirichletCondition[ψ[x]==0,x==xdagger]
};
(*Boundary conditions for the penalized PDE*)
bc={
DirichletCondition[ψ[y,x]==V0[y]-(a y-b),x==0],
DirichletCondition[ψ[y,x]==V00[x]-Phi00[x],y==0],
DirichletCondition[ψ[y,x]==0,x==xdagger],
DirichletCondition[ψ[y,x]==0,y==ydagger]};

(*DIFFERENTIAL EQUATIONS*)
(*Nonlinear ODE for the benchmark penalized problem in y*)
nonlinearODE0=L0[ψ[y],y]+(1/ϵ)Max[- ψ[y];0]==-L0[a*y-b,y];
(*Nonlinear ODE for the benchmark penalized problem in x*)
nonlinearODE00=L00[ψ[x],x]+(1/ϵ)Max[- ψ[x];0]==-L00[Phi00[x],x];
(*Nonlinear PDE for the two-dimensional penalized problem*)
nonlinearPDE=L[ψ[y,x],y,x]+(1/ϵ)Max[ -ψ[y,x];0]==-L[Phi[y,x],y,x];

(*NONLINEAR ODE/PDE SOLVERs*)
ψ0=NDSolveValue[
{nonlinearODE0,bc0},ψ,{y,0,ydagger},
Method->"FiniteElement"];
ψ00=NDSolveValue[
{nonlinearODE00,bc00},ψ,{x,0,xdagger},
Method->"FiniteElement"];
V00e[x_]:=ψ00[x]+Phi00[x];
ψ=NDSolveValue[
{nonlinearPDE,bc},ψ,{y,x}∈Omega,
Method->"FiniteElement"
];
];
Export["ψσx1.wxf",ψ];

Comparative statics graph for (Subscript[σ, x])

ψbaseline=Import["ψbaseline.wxf"];
xdagger=Import["xdagger.wxf"];
ydagger=Import["ydagger.wxf"];
ψ0=Import["ψσx0.wxf"];
ψ1=Import["ψσx1.wxf"];

Fig4b1=ContourPlot[ψ0[y,x]==(1/(10^4)),{x,xsmall,xdagger},{y,0,ydagger},
ContourStyle->{Directive[Black,Dotted,Thickness[0.005]]},Frame->{True,True,False,False},
FrameTicks->{{{1,2,3,4},None},{{0,50,100,150,200,250},None}},
PlotRange->{{0,xlarge},{0,ylarge}},
FrameLabel->{"Property price x","Fare revenues y","",""}];Fig4b2=ContourPlot[ψbaseline[y,x]==(1/(10^4)),{x,xsmall,xdagger},{y,0,ydagger},Frame->True,
ContourStyle->{Directive[Black,Plain,Thickness[0.005]]}];
Fig4b3=ContourPlot[ψ1[y,x]==(1/(10^4)),{x,xsmall,xdagger},{y,0,ydagger},ContourStyle->{Directive[Black,Dashed,Thickness[0.005]]}];
Figure4b=Legended[
Show[{Fig4b1,Fig4b2,Fig4b3},AspectRatio->aspectratio],
Placed[Framed@LineLegend[{Directive[Black,Dotted,Thickness[0.005]],Directive[Black,Plain,Thickness[0.005]],Directive[Black,Dashed,Thickness[0.005]]},{"(Subscript[σ, x])=.15","(Subscript[σ, x])=.20","(Subscript[σ, x])=.25"},LegendMarkerSize->15,Spacings->.2,LegendLabel->"Curve (Subscript[Γ, ϵ])(·) for"],{Center,Top}]];
Export["Figure4b.png",Figure4b];SystemOpen["Figure4b.png"]

Panel C -- Discount rate r

r=.100 (low)

Clear[ψ,xdagger,ydagger,ψ0,ψ1];
Module[
{
(*PRIMARY PARAMETERS*)
(*Process parameters*)
σx=.20, 
σy=.0675,
ρ=0.8,
μx=0.0725,
μy=0.049,
r=0.10,(*CHANGED PARAMETER*)
(*Profit parameters*)
I1=24.849,
c=.9,
I2=14.5,
K=4.53,
α=.51
},
(*SECONDARY PARAMETERS*)
(*Parameters for obstacle*)
a=(1/(r-μy));
b=(c/r)+I1;
β=α*I2-K;
γ=((-2 μx+((σx)^2)+2 (Sqrt[2 r ((σx)^2)+(((μx-(((σx)^2)/2)))^2)]))/(2 ((σx)^2)));
γtilde=((-2 μx+((σx)^2)-2 (Sqrt[2 r ((σx)^2)+(((μx-(((σx)^2)/2)))^2)]))/(2 ((σx)^2)));
δ=((-2 μy+((σy)^2)+2 (Sqrt[2 r ((σy)^2)+(((μy-(((σy)^2)/2)))^2)]))/(2 ((σy)^2)));
xstar=(γ/(γ-1))*(I2+(K/(1-α)));
yFB=(δ/(δ-1))*(b/a);
xFB=(γ/(γ-1))*((b+β)/α);
(*Parameters for benchmark problems*)
δe=((-2 μy+((σy)^2)+2 (Sqrt[-2 (-r-(1/ϵ)) ((σy)^2)+(((μy-(((σy)^2)/2)))^2)]))/(2 ((σy)^2)));
γe=((-2 μx+((σx)^2)+2 (Sqrt[-2 (-r-(1/ϵ)) ((σx)^2)+(((μx-(((σx)^2)/2)))^2)]))/(2 ((σx)^2)));
yFBe=((δ+ϵ*δe*(r-μy*δ)+ϵ*r*(δ-1))/(δ+ϵ*δe*(r-μy*δ)+ϵ*r*(δ-1)-1))*(b/a);
xFBe=((γ+ϵ*γe*(r-μx*γ)+ϵ*r*(γ-1))/(γ+ϵ*γe*(r-μx*γ)+ϵ*r*(γ-1)-1))*((b+β)/α);
(*Domain of PDE*)
ydagger=yFBe*daggertobarratio;
xdagger=xFBe*daggertobarratio;
Omega=Rectangle[{0,0},{ydagger,xdagger}];

(*INTERMEDIATE FUNCTIONS*)
(*Obstacles*)
Phi00[x_]:=-b + ({{\[Piecewise], ({{((((x/xstar)))^γ)*(α*xstar-β), 0<x<=xstar}, {α*x-β, x>=xstar }})}});
Phi[y_,x_]:=a*y-b + ({{\[Piecewise], ({{((((x/xstar)))^γ)*(α*xstar-β), 0<x<=xstar}, {α*x-β, x>=xstar }})}});
(*Define the function V^0(y)*)
V0[y_]:=
({{\[Piecewise], ({{b/(δ-1)*((((y/yFB)))^δ), 0<y<yFB}, {a*y-b, y>=yFB}})}});
(*Define the function V_0(x)*)
V00[x_]:=
({{\[Piecewise], ({{(α*Max[xstar,xFB]-(β+b))*((((x/(Max[xstar,xFB]))))^γ), 0<=x<Max[xstar,xFB]}, {α*x-(β+b), x>=Max[xstar,xFB]}})}});

(*DIFFERENTIAL OPERATORS*)
(*Differential operator of the ODE in y*)
L0[v_,y_]:=(1/2)((σy)^2) (y^2) D[v,{y,2}]+μy y D[v,y]-r v;
(*Differential operator of the ODE in x*)
L00[v_,x_]:=(1/2)((σx)^2) (x^2) D[v,{x,2}]+μx x D[v,x]-r v;
(*Differential operator L*)
L[v_,y_,x_]:=(1/2)((σy)^2) (y^2) D[v,{y,2}]+(1/2) ((σx)^2)(x^2) D[v,{x,2}]+ρ σy σx y x D[v,y,x]+μy y D[v,y]+μx x D[v,x]-r v;

(*BOUNDARY CONDITIONS OF THE BENCHMARK PENALIZED ODEs*)
(*Boundary conditions for the penalized ODE in y*)
bc0={
DirichletCondition[ψ[y]==b,y==0],DirichletCondition[ψ[y]==0,y==ydagger]
};
(*Boundary conditions for the penalized ODE in x*)
bc00={
DirichletCondition[ψ[x]==b,x==0],DirichletCondition[ψ[x]==0,x==xdagger]
};
(*Boundary conditions for the penalized PDE*)
bc={
DirichletCondition[ψ[y,x]==V0[y]-(a y-b),x==0],
DirichletCondition[ψ[y,x]==V00[x]-Phi00[x],y==0],
DirichletCondition[ψ[y,x]==0,x==xdagger],
DirichletCondition[ψ[y,x]==0,y==ydagger]};

(*DIFFERENTIAL EQUATIONS*)
(*Nonlinear ODE for the benchmark penalized problem in y*)
nonlinearODE0=L0[ψ[y],y]+(1/ϵ)Max[- ψ[y];0]==-L0[a*y-b,y];
(*Nonlinear ODE for the benchmark penalized problem in x*)
nonlinearODE00=L00[ψ[x],x]+(1/ϵ)Max[- ψ[x];0]==-L00[Phi00[x],x];
(*Nonlinear PDE for the two-dimensional penalized problem*)
nonlinearPDE=L[ψ[y,x],y,x]+(1/ϵ)Max[ -ψ[y,x];0]==-L[Phi[y,x],y,x];

(*NONLINEAR ODE/PDE SOLVERs*)
ψ0=NDSolveValue[
{nonlinearODE0,bc0},ψ,{y,0,ydagger},
Method->"FiniteElement"];
ψ00=NDSolveValue[
{nonlinearODE00,bc00},ψ,{x,0,xdagger},
Method->"FiniteElement"];
V00e[x_]:=ψ00[x]+Phi00[x];
ψ=NDSolveValue[
{nonlinearPDE,bc},ψ,{y,x}∈Omega,
Method->"FiniteElement"
];
];
Export["ψr0.wxf",ψ];

r=.1200 (high)

Clear[ψ,xdagger,ydagger,ψ0,ψ1];
Module[
{
(*PRIMARY PARAMETERS*)
(*Process parameters*)
σx=.20,
σy=.0675,
ρ=0.8,
μx=0.0725,
μy=0.049,
r=0.1200, (*CHANGED PARAMETER*)
(*Profit parameters*)
I1=24.849,
c=.9,
I2=14.5,
K=4.53,
α=.51
},
(*SECONDARY PARAMETERS*)
(*Parameters for obstacle*)
a=(1/(r-μy));
b=(c/r)+I1;
β=α*I2-K;
γ=((-2 μx+((σx)^2)+2 (Sqrt[2 r ((σx)^2)+(((μx-(((σx)^2)/2)))^2)]))/(2 ((σx)^2)));
γtilde=((-2 μx+((σx)^2)-2 (Sqrt[2 r ((σx)^2)+(((μx-(((σx)^2)/2)))^2)]))/(2 ((σx)^2)));
δ=((-2 μy+((σy)^2)+2 (Sqrt[2 r ((σy)^2)+(((μy-(((σy)^2)/2)))^2)]))/(2 ((σy)^2)));
xstar=(γ/(γ-1))*(I2+(K/(1-α)));
yFB=(δ/(δ-1))*(b/a);
xFB=(γ/(γ-1))*((b+β)/α);
(*Parameters for benchmark problems*)
δe=((-2 μy+((σy)^2)+2 (Sqrt[-2 (-r-(1/ϵ)) ((σy)^2)+(((μy-(((σy)^2)/2)))^2)]))/(2 ((σy)^2)));
γe=((-2 μx+((σx)^2)+2 (Sqrt[-2 (-r-(1/ϵ)) ((σx)^2)+(((μx-(((σx)^2)/2)))^2)]))/(2 ((σx)^2)));
yFBe=((δ+ϵ*δe*(r-μy*δ)+ϵ*r*(δ-1))/(δ+ϵ*δe*(r-μy*δ)+ϵ*r*(δ-1)-1))*(b/a);
xFBe=((γ+ϵ*γe*(r-μx*γ)+ϵ*r*(γ-1))/(γ+ϵ*γe*(r-μx*γ)+ϵ*r*(γ-1)-1))*((b+β)/α);
(*Domain of PDE*)
ydagger=yFBe*daggertobarratio;
xdagger=xFBe*daggertobarratio;
Omega=Rectangle[{0,0},{ydagger,xdagger}];

(*INTERMEDIATE FUNCTIONS*)
(*Obstacles*)
Phi00[x_]:=-b + ({{\[Piecewise], ({{((((x/xstar)))^γ)*(α*xstar-β), 0<x<=xstar}, {α*x-β, x>=xstar }})}});
Phi[y_,x_]:=a*y-b + ({{\[Piecewise], ({{((((x/xstar)))^γ)*(α*xstar-β), 0<x<=xstar}, {α*x-β, x>=xstar }})}});
(*Define the function V^0(y)*)
V0[y_]:=
({{\[Piecewise], ({{b/(δ-1)*((((y/yFB)))^δ), 0<y<yFB}, {a*y-b, y>=yFB}})}});
(*Define the function V_0(x)*)
V00[x_]:=
({{\[Piecewise], ({{(α*Max[xstar,xFB]-(β+b))*((((x/(Max[xstar,xFB]))))^γ), 0<=x<Max[xstar,xFB]}, {α*x-(β+b), x>=Max[xstar,xFB]}})}});

(*DIFFERENTIAL OPERATORS*)
(*Differential operator of the ODE in y*)
L0[v_,y_]:=(1/2)((σy)^2) (y^2) D[v,{y,2}]+μy y D[v,y]-r v;
(*Differential operator of the ODE in x*)
L00[v_,x_]:=(1/2)((σx)^2) (x^2) D[v,{x,2}]+μx x D[v,x]-r v;
(*Differential operator L*)
L[v_,y_,x_]:=(1/2)((σy)^2) (y^2) D[v,{y,2}]+(1/2) ((σx)^2)(x^2) D[v,{x,2}]+ρ σy σx y x D[v,y,x]+μy y D[v,y]+μx x D[v,x]-r v;

(*BOUNDARY CONDITIONS OF THE BENCHMARK PENALIZED ODEs*)
(*Boundary conditions for the penalized ODE in y*)
bc0={
DirichletCondition[ψ[y]==b,y==0],DirichletCondition[ψ[y]==0,y==ydagger]
};
(*Boundary conditions for the penalized ODE in x*)
bc00={
DirichletCondition[ψ[x]==b,x==0],DirichletCondition[ψ[x]==0,x==xdagger]
};
(*Boundary conditions for the penalized PDE*)
bc={
DirichletCondition[ψ[y,x]==V0[y]-(a y-b),x==0],
DirichletCondition[ψ[y,x]==V00[x]-Phi00[x],y==0],
DirichletCondition[ψ[y,x]==0,x==xdagger],
DirichletCondition[ψ[y,x]==0,y==ydagger]};

(*DIFFERENTIAL EQUATIONS*)
(*Nonlinear ODE for the benchmark penalized problem in y*)
nonlinearODE0=L0[ψ[y],y]+(1/ϵ)Max[- ψ[y];0]==-L0[a*y-b,y];
(*Nonlinear ODE for the benchmark penalized problem in x*)
nonlinearODE00=L00[ψ[x],x]+(1/ϵ)Max[- ψ[x];0]==-L00[Phi00[x],x];
(*Nonlinear PDE for the two-dimensional penalized problem*)
nonlinearPDE=L[ψ[y,x],y,x]+(1/ϵ)Max[ -ψ[y,x];0]==-L[Phi[y,x],y,x];

(*NONLINEAR ODE/PDE SOLVERs*)
ψ0=NDSolveValue[
{nonlinearODE0,bc0},ψ,{y,0,ydagger},
Method->"FiniteElement"];
ψ00=NDSolveValue[
{nonlinearODE00,bc00},ψ,{x,0,xdagger},
Method->"FiniteElement"];
V00e[x_]:=ψ00[x]+Phi00[x];
ψ=NDSolveValue[
{nonlinearPDE,bc},ψ,{y,x}∈Omega,
Method->"FiniteElement"
];
];
Export["ψr1.wxf",ψ];

Comparative statics graph for r

ψbaseline=Import["ψbaseline.wxf"];
xdagger=Import["xdagger.wxf"];
ydagger=Import["ydagger.wxf"];
ψ0=Import["ψr0.wxf"];
ψ1=Import["ψr1.wxf"];
Fig4c1=ContourPlot[ψ0[y,x]==(1/(10^4)),{x,xsmall,xdagger},{y,0,ydagger},
ContourStyle->{Directive[Black,Dotted,Thickness[0.005]]},Frame->{True,True,False,False},
FrameTicks->{{{1,2,3,4},None},{{0,50,100,150,200,250},None}},
PlotRange->{{0,xlarge},{0,ylarge}},
FrameLabel->{"Property price x","Fare revenues y","",""}];Fig4c2=ContourPlot[ψbaseline[y,x]==(1/(10^4)),{x,xsmall,xdagger},{y,0,ydagger},Frame->True,
ContourStyle->{Directive[Black,Plain,Thickness[0.005]]}];
Fig4c3=ContourPlot[ψ1[y,x]==(1/(10^4)),{x,xsmall*2,xdagger},{y,0,ydagger},ContourStyle->{Directive[Black,Dashed,Thickness[0.005]]}];
Figure4c=Legended[
Show[{Fig4c1,Fig4c2,Fig4c3},AspectRatio->aspectratio],Placed[Framed@LineLegend[{Directive[Black,Dotted,Thickness[0.005]],Directive[Black,Plain,Thickness[0.005]],Directive[Black,Dashed,Thickness[0.005]]},{"r=.1000","r=.1133","r=.1200"},LegendMarkerSize->15,Spacings->.2,LegendLabel->"Curve (Subscript[Γ, ϵ])(·) for"],{Center,Top}]];
Export["Figure4c.png",Figure4c];SystemOpen["Figure4c.png"]

Panel D -- Growth rate (Subscript[μ, y])

(Subscript[μ, y])=.020 (low)

Clear[ψ,xdagger,ydagger,ψ0,ψ1];
Module[
{
(*PRIMARY PARAMETERS*)
(*Process parameters*)
σx=.2,
σy=.0675,
ρ=0.8,
μx=.0725, 
μy=.020,(*CHANGED PARAMETER*)
r=0.1133,
(*Profit parameters*)
I1=24.849,
c=.9,
I2=14.5,
K=4.53,
α=.51
},
(*SECONDARY PARAMETERS*)
(*Parameters for obstacle*)
a=(1/(r-μy));
b=(c/r)+I1;
β=α*I2-K;
γ=((-2 μx+((σx)^2)+2 (Sqrt[2 r ((σx)^2)+(((μx-(((σx)^2)/2)))^2)]))/(2 ((σx)^2)));
γtilde=((-2 μx+((σx)^2)-2 (Sqrt[2 r ((σx)^2)+(((μx-(((σx)^2)/2)))^2)]))/(2 ((σx)^2)));
δ=((-2 μy+((σy)^2)+2 (Sqrt[2 r ((σy)^2)+(((μy-(((σy)^2)/2)))^2)]))/(2 ((σy)^2)));
xstar=(γ/(γ-1))*(I2+(K/(1-α)));
yFB=(δ/(δ-1))*(b/a);
xFB=(γ/(γ-1))*((b+β)/α);
(*Parameters for benchmark problems*)
δe=((-2 μy+((σy)^2)+2 (Sqrt[-2 (-r-(1/ϵ)) ((σy)^2)+(((μy-(((σy)^2)/2)))^2)]))/(2 ((σy)^2)));
γe=((-2 μx+((σx)^2)+2 (Sqrt[-2 (-r-(1/ϵ)) ((σx)^2)+(((μx-(((σx)^2)/2)))^2)]))/(2 ((σx)^2)));
yFBe=((δ+ϵ*δe*(r-μy*δ)+ϵ*r*(δ-1))/(δ+ϵ*δe*(r-μy*δ)+ϵ*r*(δ-1)-1))*(b/a);
xFBe=((γ+ϵ*γe*(r-μx*γ)+ϵ*r*(γ-1))/(γ+ϵ*γe*(r-μx*γ)+ϵ*r*(γ-1)-1))*((b+β)/α);
(*Domain of PDE*)
ydagger=yFBe*daggertobarratio;
xdagger=xFBe*daggertobarratio;
Omega=Rectangle[{0,0},{ydagger,xdagger}];

(*INTERMEDIATE FUNCTIONS*)
(*Obstacles*)
Phi00[x_]:=-b + ({{\[Piecewise], ({{((((x/xstar)))^γ)*(α*xstar-β), 0<x<=xstar}, {α*x-β, x>=xstar }})}});
Phi[y_,x_]:=a*y-b + ({{\[Piecewise], ({{((((x/xstar)))^γ)*(α*xstar-β), 0<x<=xstar}, {α*x-β, x>=xstar }})}});
(*Define the function V^0(y)*)
V0[y_]:=
({{\[Piecewise], ({{b/(δ-1)*((((y/yFB)))^δ), 0<y<yFB}, {a*y-b, y>=yFB}})}});
(*Define the function V_0(x)*)
V00[x_]:=
({{\[Piecewise], ({{(α*Max[xstar,xFB]-(β+b))*((((x/(Max[xstar,xFB]))))^γ), 0<=x<Max[xstar,xFB]}, {α*x-(β+b), x>=Max[xstar,xFB]}})}});

(*DIFFERENTIAL OPERATORS*)
(*Differential operator of the ODE in y*)
L0[v_,y_]:=(1/2)((σy)^2) (y^2) D[v,{y,2}]+μy y D[v,y]-r v;
(*Differential operator of the ODE in x*)
L00[v_,x_]:=(1/2)((σx)^2) (x^2) D[v,{x,2}]+μx x D[v,x]-r v;
(*Differential operator L*)
L[v_,y_,x_]:=(1/2)((σy)^2) (y^2) D[v,{y,2}]+(1/2) ((σx)^2)(x^2) D[v,{x,2}]+ρ σy σx y x D[v,y,x]+μy y D[v,y]+μx x D[v,x]-r v;

(*BOUNDARY CONDITIONS OF THE BENCHMARK PENALIZED ODEs*)
(*Boundary conditions for the penalized ODE in y*)
bc0={
DirichletCondition[ψ[y]==b,y==0],DirichletCondition[ψ[y]==0,y==ydagger]
};
(*Boundary conditions for the penalized ODE in x*)
bc00={
DirichletCondition[ψ[x]==b,x==0],DirichletCondition[ψ[x]==0,x==xdagger]
};
(*Boundary conditions for the penalized PDE*)
bc={
DirichletCondition[ψ[y,x]==V0[y]-(a y-b),x==0],
DirichletCondition[ψ[y,x]==V00[x]-Phi00[x],y==0],
DirichletCondition[ψ[y,x]==0,x==xdagger],
DirichletCondition[ψ[y,x]==0,y==ydagger]};

(*DIFFERENTIAL EQUATIONS*)
(*Nonlinear ODE for the benchmark penalized problem in y*)
nonlinearODE0=L0[ψ[y],y]+(1/ϵ)Max[- ψ[y];0]==-L0[a*y-b,y];
(*Nonlinear ODE for the benchmark penalized problem in x*)
nonlinearODE00=L00[ψ[x],x]+(1/ϵ)Max[- ψ[x];0]==-L00[Phi00[x],x];
(*Nonlinear PDE for the two-dimensional penalized problem*)
nonlinearPDE=L[ψ[y,x],y,x]+(1/ϵ)Max[ -ψ[y,x];0]==-L[Phi[y,x],y,x];

(*NONLINEAR ODE/PDE SOLVERs*)
ψ0=NDSolveValue[
{nonlinearODE0,bc0},ψ,{y,0,ydagger},
Method->"FiniteElement"];
ψ00=NDSolveValue[
{nonlinearODE00,bc00},ψ,{x,0,xdagger},
Method->"FiniteElement"];
V00e[x_]:=ψ00[x]+Phi00[x];
ψ=NDSolveValue[
{nonlinearPDE,bc},ψ,{y,x}∈Omega,
Method->"FiniteElement"
];
];
Export["ψμy0.wxf",ψ];

(Subscript[μ, y])=.0600 (high)

Clear[ψ,xdagger,ydagger];
Module[
{
(*PRIMARY PARAMETERS*)
(*Process parameters*)
σx=.2,
σy=.0675,
ρ=0.8,
μx=.0725, 
μy=.060,(*CHANGED PARAMETER*)
r=0.1133,
(*Profit parameters*)
I1=24.849,
c=.9,
I2=14.5,
K=4.53,
α=.51
},
(*SECONDARY PARAMETERS*)
(*Parameters for obstacle*)
a=(1/(r-μy));
b=(c/r)+I1;
β=α*I2-K;
γ=((-2 μx+((σx)^2)+2 (Sqrt[2 r ((σx)^2)+(((μx-(((σx)^2)/2)))^2)]))/(2 ((σx)^2)));
γtilde=((-2 μx+((σx)^2)-2 (Sqrt[2 r ((σx)^2)+(((μx-(((σx)^2)/2)))^2)]))/(2 ((σx)^2)));
δ=((-2 μy+((σy)^2)+2 (Sqrt[2 r ((σy)^2)+(((μy-(((σy)^2)/2)))^2)]))/(2 ((σy)^2)));
xstar=(γ/(γ-1))*(I2+(K/(1-α)));
yFB=(δ/(δ-1))*(b/a);
xFB=(γ/(γ-1))*((b+β)/α);
(*Parameters for benchmark problems*)
δe=((-2 μy+((σy)^2)+2 (Sqrt[-2 (-r-(1/ϵ)) ((σy)^2)+(((μy-(((σy)^2)/2)))^2)]))/(2 ((σy)^2)));
γe=((-2 μx+((σx)^2)+2 (Sqrt[-2 (-r-(1/ϵ)) ((σx)^2)+(((μx-(((σx)^2)/2)))^2)]))/(2 ((σx)^2)));
yFBe=((δ+ϵ*δe*(r-μy*δ)+ϵ*r*(δ-1))/(δ+ϵ*δe*(r-μy*δ)+ϵ*r*(δ-1)-1))*(b/a);
xFBe=((γ+ϵ*γe*(r-μx*γ)+ϵ*r*(γ-1))/(γ+ϵ*γe*(r-μx*γ)+ϵ*r*(γ-1)-1))*((b+β)/α);
(*Domain of PDE*)
ydagger=yFBe*daggertobarratio;
xdagger=xFBe*daggertobarratio;
Omega=Rectangle[{0,0},{ydagger,xdagger}];

(*INTERMEDIATE FUNCTIONS*)
(*Obstacles*)
Phi00[x_]:=-b + ({{\[Piecewise], ({{((((x/xstar)))^γ)*(α*xstar-β), 0<x<=xstar}, {α*x-β, x>=xstar }})}});
Phi[y_,x_]:=a*y-b + ({{\[Piecewise], ({{((((x/xstar)))^γ)*(α*xstar-β), 0<x<=xstar}, {α*x-β, x>=xstar }})}});
(*Define the function V^0(y)*)
V0[y_]:=
({{\[Piecewise], ({{b/(δ-1)*((((y/yFB)))^δ), 0<y<yFB}, {a*y-b, y>=yFB}})}});
(*Define the function V_0(x)*)
V00[x_]:=
({{\[Piecewise], ({{(α*Max[xstar,xFB]-(β+b))*((((x/(Max[xstar,xFB]))))^γ), 0<=x<Max[xstar,xFB]}, {α*x-(β+b), x>=Max[xstar,xFB]}})}});

(*DIFFERENTIAL OPERATORS*)
(*Differential operator of the ODE in y*)
L0[v_,y_]:=(1/2)((σy)^2) (y^2) D[v,{y,2}]+μy y D[v,y]-r v;
(*Differential operator of the ODE in x*)
L00[v_,x_]:=(1/2)((σx)^2) (x^2) D[v,{x,2}]+μx x D[v,x]-r v;
(*Differential operator L*)
L[v_,y_,x_]:=(1/2)((σy)^2) (y^2) D[v,{y,2}]+(1/2) ((σx)^2)(x^2) D[v,{x,2}]+ρ σy σx y x D[v,y,x]+μy y D[v,y]+μx x D[v,x]-r v;

(*BOUNDARY CONDITIONS OF THE BENCHMARK PENALIZED ODEs*)
(*Boundary conditions for the penalized ODE in y*)
bc0={
DirichletCondition[ψ[y]==b,y==0],DirichletCondition[ψ[y]==0,y==ydagger]
};
(*Boundary conditions for the penalized ODE in x*)
bc00={
DirichletCondition[ψ[x]==b,x==0],DirichletCondition[ψ[x]==0,x==xdagger]
};
(*Boundary conditions for the penalized PDE*)
bc={
DirichletCondition[ψ[y,x]==V0[y]-(a y-b),x==0],
DirichletCondition[ψ[y,x]==V00[x]-Phi00[x],y==0],
DirichletCondition[ψ[y,x]==0,x==xdagger],
DirichletCondition[ψ[y,x]==0,y==ydagger]};

(*DIFFERENTIAL EQUATIONS*)
(*Nonlinear ODE for the benchmark penalized problem in y*)
nonlinearODE0=L0[ψ[y],y]+(1/ϵ)Max[- ψ[y];0]==-L0[a*y-b,y];
(*Nonlinear ODE for the benchmark penalized problem in x*)
nonlinearODE00=L00[ψ[x],x]+(1/ϵ)Max[- ψ[x];0]==-L00[Phi00[x],x];
(*Nonlinear PDE for the two-dimensional penalized problem*)
nonlinearPDE=L[ψ[y,x],y,x]+(1/ϵ)Max[ -ψ[y,x];0]==-L[Phi[y,x],y,x];

(*NONLINEAR ODE/PDE SOLVERs*)
ψ0=NDSolveValue[
{nonlinearODE0,bc0},ψ,{y,0,ydagger},
Method->"FiniteElement"];
ψ00=NDSolveValue[
{nonlinearODE00,bc00},ψ,{x,0,xdagger},
Method->"FiniteElement"];
V00e[x_]:=ψ00[x]+Phi00[x];
ψ=NDSolveValue[
{nonlinearPDE,bc},ψ,{y,x}∈Omega,
Method->"FiniteElement"
];
];
Export["ψμy1.wxf",ψ];

Comparative statics graph for (Subscript[μ, y])

ψbaseline=Import["ψbaseline.wxf"];
xdagger=Import["xdagger.wxf"];
ydagger=Import["ydagger.wxf"];
ψ0=Import["ψμy0.wxf"];
ψ1=Import["ψμy1.wxf"];
Fig4d1=ContourPlot[ψ0[y,x]==(1/(10^4)),{x,xsmall,xdagger},{y,0,ydagger},
ContourStyle->{Directive[Black,Dotted,Thickness[0.005]]},
Frame->{True,True,False,False},
PlotRange->{{0,xlarge},{0,ylarge}},
FrameTicks->{{{1,2,3,4},None},{{0,50,100,150,200,250},None}},FrameLabel->{"Property price x","Fare revenues y","",""}
];
Fig4d2=ContourPlot[ψbaseline[y,x]==(1/(10^4)),{x,xsmall,xdagger},{y,0,ydagger},Frame->True,
ContourStyle->{Directive[Black,Plain,Thickness[0.005]]}];
Fig4d3=ContourPlot[ψ1[y,x]==(1/(10^4)),{x,xsmall,xdagger},{y,0,ydagger},ContourStyle->{Directive[Black,Dashed,Thickness[0.005]]}];
Figure4d=Legended[
Show[{Fig4d1,Fig4d2,Fig4d3},AspectRatio->aspectratio],Placed[Framed@LineLegend[{Directive[Black,Dotted,Thickness[0.005]],Directive[Black,Plain,Thickness[0.005]],Directive[Black,Dashed,Thickness[0.005]]},{"(Subscript[μ, y])=.020","(Subscript[μ, y])=.049","(Subscript[μ, y])=.060"},LegendMarkerSize->15,Spacings->.2,LegendLabel->"Curve (Subscript[Γ, ϵ])(·) for"],{Center,Top}]];
Export["Figure4d.png",Figure4d];SystemOpen["Figure4d.png"]

Panel E -- Volatility (Subscript[σ, y])

(Subscript[σ, y])=.0300 (low)

Clear[ψ,xdagger,ydagger,ψ0,ψ1];
Module[
{
(*PRIMARY PARAMETERS*)
(*Process parameters*)
σx=.20,
σy=.0300, (*CHANGED PARAMETER*)
ρ=0.8,
μx=0.0725,
μy=0.049,
r=0.1133,
(*Profit parameters*)
I1=24.849,
c=.9,
I2=14.5,
K=4.53,
α=.51
},
(*SECONDARY PARAMETERS*)
(*Parameters for obstacle*)
a=(1/(r-μy));
b=(c/r)+I1;
β=α*I2-K;
γ=((-2 μx+((σx)^2)+2 (Sqrt[2 r ((σx)^2)+(((μx-(((σx)^2)/2)))^2)]))/(2 ((σx)^2)));
γtilde=((-2 μx+((σx)^2)-2 (Sqrt[2 r ((σx)^2)+(((μx-(((σx)^2)/2)))^2)]))/(2 ((σx)^2)));
δ=((-2 μy+((σy)^2)+2 (Sqrt[2 r ((σy)^2)+(((μy-(((σy)^2)/2)))^2)]))/(2 ((σy)^2)));
xstar=(γ/(γ-1))*(I2+(K/(1-α)));
yFB=(δ/(δ-1))*(b/a);
xFB=(γ/(γ-1))*((b+β)/α);
(*Parameters for benchmark problems*)
δe=((-2 μy+((σy)^2)+2 (Sqrt[-2 (-r-(1/ϵ)) ((σy)^2)+(((μy-(((σy)^2)/2)))^2)]))/(2 ((σy)^2)));
γe=((-2 μx+((σx)^2)+2 (Sqrt[-2 (-r-(1/ϵ)) ((σx)^2)+(((μx-(((σx)^2)/2)))^2)]))/(2 ((σx)^2)));
yFBe=((δ+ϵ*δe*(r-μy*δ)+ϵ*r*(δ-1))/(δ+ϵ*δe*(r-μy*δ)+ϵ*r*(δ-1)-1))*(b/a);
xFBe=((γ+ϵ*γe*(r-μx*γ)+ϵ*r*(γ-1))/(γ+ϵ*γe*(r-μx*γ)+ϵ*r*(γ-1)-1))*((b+β)/α);
(*Domain of PDE*)
ydagger=yFBe*daggertobarratio;
xdagger=xFBe*daggertobarratio;
Omega=Rectangle[{0,0},{ydagger,xdagger}];

(*INTERMEDIATE FUNCTIONS*)
(*Obstacles*)
Phi00[x_]:=-b + ({{\[Piecewise], ({{((((x/xstar)))^γ)*(α*xstar-β), 0<x<=xstar}, {α*x-β, x>=xstar }})}});
Phi[y_,x_]:=a*y-b + ({{\[Piecewise], ({{((((x/xstar)))^γ)*(α*xstar-β), 0<x<=xstar}, {α*x-β, x>=xstar }})}});
(*Define the function V^0(y)*)
V0[y_]:=
({{\[Piecewise], ({{b/(δ-1)*((((y/yFB)))^δ), 0<y<yFB}, {a*y-b, y>=yFB}})}});
(*Define the function V_0(x)*)
V00[x_]:=
({{\[Piecewise], ({{(α*Max[xstar,xFB]-(β+b))*((((x/(Max[xstar,xFB]))))^γ), 0<=x<Max[xstar,xFB]}, {α*x-(β+b), x>=Max[xstar,xFB]}})}});

(*DIFFERENTIAL OPERATORS*)
(*Differential operator of the ODE in y*)
L0[v_,y_]:=(1/2)((σy)^2) (y^2) D[v,{y,2}]+μy y D[v,y]-r v;
(*Differential operator of the ODE in x*)
L00[v_,x_]:=(1/2)((σx)^2) (x^2) D[v,{x,2}]+μx x D[v,x]-r v;
(*Differential operator L*)
L[v_,y_,x_]:=(1/2)((σy)^2) (y^2) D[v,{y,2}]+(1/2) ((σx)^2)(x^2) D[v,{x,2}]+ρ σy σx y x D[v,y,x]+μy y D[v,y]+μx x D[v,x]-r v;

(*BOUNDARY CONDITIONS OF THE BENCHMARK PENALIZED ODEs*)
(*Boundary conditions for the penalized ODE in y*)
bc0={
DirichletCondition[ψ[y]==b,y==0],DirichletCondition[ψ[y]==0,y==ydagger]
};
(*Boundary conditions for the penalized ODE in x*)
bc00={
DirichletCondition[ψ[x]==b,x==0],DirichletCondition[ψ[x]==0,x==xdagger]
};
(*Boundary conditions for the penalized PDE*)
bc={
DirichletCondition[ψ[y,x]==V0[y]-(a y-b),x==0],
DirichletCondition[ψ[y,x]==V00[x]-Phi00[x],y==0],
DirichletCondition[ψ[y,x]==0,x==xdagger],
DirichletCondition[ψ[y,x]==0,y==ydagger]};

(*DIFFERENTIAL EQUATIONS*)
(*Nonlinear ODE for the benchmark penalized problem in y*)
nonlinearODE0=L0[ψ[y],y]+(1/ϵ)Max[- ψ[y];0]==-L0[a*y-b,y];
(*Nonlinear ODE for the benchmark penalized problem in x*)
nonlinearODE00=L00[ψ[x],x]+(1/ϵ)Max[- ψ[x];0]==-L00[Phi00[x],x];
(*Nonlinear PDE for the two-dimensional penalized problem*)
nonlinearPDE=L[ψ[y,x],y,x]+(1/ϵ)Max[ -ψ[y,x];0]==-L[Phi[y,x],y,x];

(*NONLINEAR ODE/PDE SOLVERs*)
ψ0=NDSolveValue[
{nonlinearODE0,bc0},ψ,{y,0,ydagger},
Method->"FiniteElement"];
ψ00=NDSolveValue[
{nonlinearODE00,bc00},ψ,{x,0,xdagger},
Method->"FiniteElement"];
V00e[x_]:=ψ00[x]+Phi00[x];
ψ=NDSolveValue[
{nonlinearPDE,bc},ψ,{y,x}∈Omega,
Method->"FiniteElement"
];
];
Export["ψσy0.wxf",ψ];

(Subscript[σ, y])=.1000 (high)

Clear[ψ,xdagger,ydagger];
Module[
{
(*PRIMARY PARAMETERS*)
(*Process parameters*)
σx=.20,
σy=.1000, (*CHANGED PARAMETER*)
ρ=0.8,
μx=0.0725,
μy=0.049,
r=0.1133,
(*Profit parameters*)
I1=24.849,
c=.9,
I2=14.5,
K=4.53,
α=.51
},
(*SECONDARY PARAMETERS*)
(*Parameters for obstacle*)
a=(1/(r-μy));
b=(c/r)+I1;
β=α*I2-K;
γ=((-2 μx+((σx)^2)+2 (Sqrt[2 r ((σx)^2)+(((μx-(((σx)^2)/2)))^2)]))/(2 ((σx)^2)));
γtilde=((-2 μx+((σx)^2)-2 (Sqrt[2 r ((σx)^2)+(((μx-(((σx)^2)/2)))^2)]))/(2 ((σx)^2)));
δ=((-2 μy+((σy)^2)+2 (Sqrt[2 r ((σy)^2)+(((μy-(((σy)^2)/2)))^2)]))/(2 ((σy)^2)));
xstar=(γ/(γ-1))*(I2+(K/(1-α)));
yFB=(δ/(δ-1))*(b/a);
xFB=(γ/(γ-1))*((b+β)/α);
(*Parameters for benchmark problems*)
δe=((-2 μy+((σy)^2)+2 (Sqrt[-2 (-r-(1/ϵ)) ((σy)^2)+(((μy-(((σy)^2)/2)))^2)]))/(2 ((σy)^2)));
γe=((-2 μx+((σx)^2)+2 (Sqrt[-2 (-r-(1/ϵ)) ((σx)^2)+(((μx-(((σx)^2)/2)))^2)]))/(2 ((σx)^2)));
yFBe=((δ+ϵ*δe*(r-μy*δ)+ϵ*r*(δ-1))/(δ+ϵ*δe*(r-μy*δ)+ϵ*r*(δ-1)-1))*(b/a);
xFBe=((γ+ϵ*γe*(r-μx*γ)+ϵ*r*(γ-1))/(γ+ϵ*γe*(r-μx*γ)+ϵ*r*(γ-1)-1))*((b+β)/α);
(*Domain of PDE*)
ydagger=yFBe*daggertobarratio;
xdagger=xFBe*daggertobarratio;
Omega=Rectangle[{0,0},{ydagger,xdagger}];

(*INTERMEDIATE FUNCTIONS*)
(*Obstacles*)
Phi00[x_]:=-b + ({{\[Piecewise], ({{((((x/xstar)))^γ)*(α*xstar-β), 0<x<=xstar}, {α*x-β, x>=xstar }})}});
Phi[y_,x_]:=a*y-b + ({{\[Piecewise], ({{((((x/xstar)))^γ)*(α*xstar-β), 0<x<=xstar}, {α*x-β, x>=xstar }})}});
(*Define the function V^0(y)*)
V0[y_]:=
({{\[Piecewise], ({{b/(δ-1)*((((y/yFB)))^δ), 0<y<yFB}, {a*y-b, y>=yFB}})}});
(*Define the function V_0(x)*)
V00[x_]:=
({{\[Piecewise], ({{(α*Max[xstar,xFB]-(β+b))*((((x/(Max[xstar,xFB]))))^γ), 0<=x<Max[xstar,xFB]}, {α*x-(β+b), x>=Max[xstar,xFB]}})}});

(*DIFFERENTIAL OPERATORS*)
(*Differential operator of the ODE in y*)
L0[v_,y_]:=(1/2)((σy)^2) (y^2) D[v,{y,2}]+μy y D[v,y]-r v;
(*Differential operator of the ODE in x*)
L00[v_,x_]:=(1/2)((σx)^2) (x^2) D[v,{x,2}]+μx x D[v,x]-r v;
(*Differential operator L*)
L[v_,y_,x_]:=(1/2)((σy)^2) (y^2) D[v,{y,2}]+(1/2) ((σx)^2)(x^2) D[v,{x,2}]+ρ σy σx y x D[v,y,x]+μy y D[v,y]+μx x D[v,x]-r v;

(*BOUNDARY CONDITIONS OF THE BENCHMARK PENALIZED ODEs*)
(*Boundary conditions for the penalized ODE in y*)
bc0={
DirichletCondition[ψ[y]==b,y==0],DirichletCondition[ψ[y]==0,y==ydagger]
};
(*Boundary conditions for the penalized ODE in x*)
bc00={
DirichletCondition[ψ[x]==b,x==0],DirichletCondition[ψ[x]==0,x==xdagger]
};
(*Boundary conditions for the penalized PDE*)
bc={
DirichletCondition[ψ[y,x]==V0[y]-(a y-b),x==0],
DirichletCondition[ψ[y,x]==V00[x]-Phi00[x],y==0],
DirichletCondition[ψ[y,x]==0,x==xdagger],
DirichletCondition[ψ[y,x]==0,y==ydagger]};

(*DIFFERENTIAL EQUATIONS*)
(*Nonlinear ODE for the benchmark penalized problem in y*)
nonlinearODE0=L0[ψ[y],y]+(1/ϵ)Max[- ψ[y];0]==-L0[a*y-b,y];
(*Nonlinear ODE for the benchmark penalized problem in x*)
nonlinearODE00=L00[ψ[x],x]+(1/ϵ)Max[- ψ[x];0]==-L00[Phi00[x],x];
(*Nonlinear PDE for the two-dimensional penalized problem*)
nonlinearPDE=L[ψ[y,x],y,x]+(1/ϵ)Max[ -ψ[y,x];0]==-L[Phi[y,x],y,x];

(*NONLINEAR ODE/PDE SOLVERs*)
ψ0=NDSolveValue[
{nonlinearODE0,bc0},ψ,{y,0,ydagger},
Method->"FiniteElement"];
ψ00=NDSolveValue[
{nonlinearODE00,bc00},ψ,{x,0,xdagger},
Method->"FiniteElement"];
V00e[x_]:=ψ00[x]+Phi00[x];
ψ=NDSolveValue[
{nonlinearPDE,bc},ψ,{y,x}∈Omega,
Method->"FiniteElement"
];
];
Export["ψσy1.wxf",ψ];

Comparative statics graph for (Subscript[σ, y])

ψbaseline=Import["ψbaseline.wxf"];
xdagger=Import["xdagger.wxf"];
ydagger=Import["ydagger.wxf"];
ψ0=Import["ψσy0.wxf"];
ψ1=Import["ψσy1.wxf"];

Fig4e1=ContourPlot[ψ0[y,x]==(1/(10^4)),{x,xsmall,xdagger},{y,0,ydagger},
ContourStyle->{Directive[Black,Dotted,Thickness[0.005]]},Frame->{True,True,False,False},
FrameTicks->{{{1,2,3,4},None},{{0,50,100,150,200,250},None}},
PlotRange->{{0,xlarge},{0,ylarge}},
FrameLabel->{"Property price x","Fare revenues y","",""}];Fig4e2=ContourPlot[ψbaseline[y,x]==(1/(10^4)),{x,xsmall,xdagger},{y,0,ydagger},Frame->True,
ContourStyle->{Directive[Black,Plain,Thickness[0.005]]}];
Fig4e3=ContourPlot[ψ1[y,x]==(1/(10^4)),{x,xsmall,xdagger},{y,0,ydagger},ContourStyle->{Directive[Black,Dashed,Thickness[0.005]]}];
Figure4e=Legended[
Show[{Fig4e1,Fig4e2,Fig4e3},AspectRatio->aspectratio],Placed[Framed@LineLegend[{Directive[Black,Dotted,Thickness[0.005]],Directive[Black,Plain,Thickness[0.005]],Directive[Black,Dashed,Thickness[0.005]]},{"(Subscript[σ, y])=.0300","(Subscript[σ, y])=.0675","(Subscript[σ, y])=.1000"},LegendMarkerSize->15,Spacings->.2,LegendLabel->"Curve (Subscript[Γ, ϵ])(·) for"],{Center,Top}]];
Export["Figure4e.png",Figure4e];SystemOpen["Figure4e.png"]

Panel F -- Correlation ρ

ρ=-.8 (negative)

Clear[ψ,xdagger,ydagger,ψ0,ψ1];
Module[
{
(*PRIMARY PARAMETERS*)
(*Process parameters*)
σx=.20,
σy=.0675,
ρ=-0.8, (*CHANGED PARAMETER*)
μx=0.0725,
μy=0.049,
r=0.1133,
(*Profit parameters*)
I1=24.849,
c=.9,
I2=14.5,
K=4.53,
α=.51
},
(*SECONDARY PARAMETERS*)
(*Parameters for obstacle*)
a=(1/(r-μy));
b=(c/r)+I1;
β=α*I2-K;
γ=((-2 μx+((σx)^2)+2 (Sqrt[2 r ((σx)^2)+(((μx-(((σx)^2)/2)))^2)]))/(2 ((σx)^2)));
γtilde=((-2 μx+((σx)^2)-2 (Sqrt[2 r ((σx)^2)+(((μx-(((σx)^2)/2)))^2)]))/(2 ((σx)^2)));
δ=((-2 μy+((σy)^2)+2 (Sqrt[2 r ((σy)^2)+(((μy-(((σy)^2)/2)))^2)]))/(2 ((σy)^2)));
xstar=(γ/(γ-1))*(I2+(K/(1-α)));
yFB=(δ/(δ-1))*(b/a);
xFB=(γ/(γ-1))*((b+β)/α);
(*Parameters for benchmark problems*)
δe=((-2 μy+((σy)^2)+2 (Sqrt[-2 (-r-(1/ϵ)) ((σy)^2)+(((μy-(((σy)^2)/2)))^2)]))/(2 ((σy)^2)));
γe=((-2 μx+((σx)^2)+2 (Sqrt[-2 (-r-(1/ϵ)) ((σx)^2)+(((μx-(((σx)^2)/2)))^2)]))/(2 ((σx)^2)));
yFBe=((δ+ϵ*δe*(r-μy*δ)+ϵ*r*(δ-1))/(δ+ϵ*δe*(r-μy*δ)+ϵ*r*(δ-1)-1))*(b/a);
xFBe=((γ+ϵ*γe*(r-μx*γ)+ϵ*r*(γ-1))/(γ+ϵ*γe*(r-μx*γ)+ϵ*r*(γ-1)-1))*((b+β)/α);
(*Domain of PDE*)
ydagger=yFBe*daggertobarratio;
xdagger=xFBe*daggertobarratio;
Omega=Rectangle[{0,0},{ydagger,xdagger}];

(*INTERMEDIATE FUNCTIONS*)
(*Obstacles*)
Phi00[x_]:=-b + ({{\[Piecewise], ({{((((x/xstar)))^γ)*(α*xstar-β), 0<x<=xstar}, {α*x-β, x>=xstar }})}});
Phi[y_,x_]:=a*y-b + ({{\[Piecewise], ({{((((x/xstar)))^γ)*(α*xstar-β), 0<x<=xstar}, {α*x-β, x>=xstar }})}});
(*Define the function V^0(y)*)
V0[y_]:=
({{\[Piecewise], ({{b/(δ-1)*((((y/yFB)))^δ), 0<y<yFB}, {a*y-b, y>=yFB}})}});
(*Define the function V_0(x)*)
V00[x_]:=
({{\[Piecewise], ({{(α*Max[xstar,xFB]-(β+b))*((((x/(Max[xstar,xFB]))))^γ), 0<=x<Max[xstar,xFB]}, {α*x-(β+b), x>=Max[xstar,xFB]}})}});

(*DIFFERENTIAL OPERATORS*)
(*Differential operator of the ODE in y*)
L0[v_,y_]:=(1/2)((σy)^2) (y^2) D[v,{y,2}]+μy y D[v,y]-r v;
(*Differential operator of the ODE in x*)
L00[v_,x_]:=(1/2)((σx)^2) (x^2) D[v,{x,2}]+μx x D[v,x]-r v;
(*Differential operator L*)
L[v_,y_,x_]:=(1/2)((σy)^2) (y^2) D[v,{y,2}]+(1/2) ((σx)^2)(x^2) D[v,{x,2}]+ρ σy σx y x D[v,y,x]+μy y D[v,y]+μx x D[v,x]-r v;

(*BOUNDARY CONDITIONS OF THE BENCHMARK PENALIZED ODEs*)
(*Boundary conditions for the penalized ODE in y*)
bc0={
DirichletCondition[ψ[y]==b,y==0],DirichletCondition[ψ[y]==0,y==ydagger]
};
(*Boundary conditions for the penalized ODE in x*)
bc00={
DirichletCondition[ψ[x]==b,x==0],DirichletCondition[ψ[x]==0,x==xdagger]
};
(*Boundary conditions for the penalized PDE*)
bc={
DirichletCondition[ψ[y,x]==V0[y]-(a y-b),x==0],
DirichletCondition[ψ[y,x]==V00[x]-Phi00[x],y==0],
DirichletCondition[ψ[y,x]==0,x==xdagger],
DirichletCondition[ψ[y,x]==0,y==ydagger]};

(*DIFFERENTIAL EQUATIONS*)
(*Nonlinear ODE for the benchmark penalized problem in y*)
nonlinearODE0=L0[ψ[y],y]+(1/ϵ)Max[- ψ[y];0]==-L0[a*y-b,y];
(*Nonlinear ODE for the benchmark penalized problem in x*)
nonlinearODE00=L00[ψ[x],x]+(1/ϵ)Max[- ψ[x];0]==-L00[Phi00[x],x];
(*Nonlinear PDE for the two-dimensional penalized problem*)
nonlinearPDE=L[ψ[y,x],y,x]+(1/ϵ)Max[ -ψ[y,x];0]==-L[Phi[y,x],y,x];

(*NONLINEAR ODE/PDE SOLVERs*)
ψ0=NDSolveValue[
{nonlinearODE0,bc0},ψ,{y,0,ydagger},
Method->"FiniteElement"];
ψ00=NDSolveValue[
{nonlinearODE00,bc00},ψ,{x,0,xdagger},
Method->"FiniteElement"];
V00e[x_]:=ψ00[x]+Phi00[x];
ψ=NDSolveValue[
{nonlinearPDE,bc},ψ,{y,x}∈Omega,
Method->"FiniteElement"
];
];
Export["ψρ0.wxf",ψ];

ρ=0 (no correlation)

Clear[ψ,xdagger,ydagger,ψ0,ψ1];
Module[
{
(*PRIMARY PARAMETERS*)
(*Process parameters*)
σx=.20,
σy=.0675,
ρ=0.0, (*CHANGED PARAMETER*)
μx=0.0725,
μy=0.049,
r=0.1133,
(*Profit parameters*)
I1=24.849,
c=.9,
I2=14.5,
K=4.53,
α=.51
},
(*SECONDARY PARAMETERS*)
(*Parameters for obstacle*)
a=(1/(r-μy));
b=(c/r)+I1;
β=α*I2-K;
γ=((-2 μx+((σx)^2)+2 (Sqrt[2 r ((σx)^2)+(((μx-(((σx)^2)/2)))^2)]))/(2 ((σx)^2)));
γtilde=((-2 μx+((σx)^2)-2 (Sqrt[2 r ((σx)^2)+(((μx-(((σx)^2)/2)))^2)]))/(2 ((σx)^2)));
δ=((-2 μy+((σy)^2)+2 (Sqrt[2 r ((σy)^2)+(((μy-(((σy)^2)/2)))^2)]))/(2 ((σy)^2)));
xstar=(γ/(γ-1))*(I2+(K/(1-α)));
yFB=(δ/(δ-1))*(b/a);
xFB=(γ/(γ-1))*((b+β)/α);
(*Parameters for benchmark problems*)
δe=((-2 μy+((σy)^2)+2 (Sqrt[-2 (-r-(1/ϵ)) ((σy)^2)+(((μy-(((σy)^2)/2)))^2)]))/(2 ((σy)^2)));
γe=((-2 μx+((σx)^2)+2 (Sqrt[-2 (-r-(1/ϵ)) ((σx)^2)+(((μx-(((σx)^2)/2)))^2)]))/(2 ((σx)^2)));
yFBe=((δ+ϵ*δe*(r-μy*δ)+ϵ*r*(δ-1))/(δ+ϵ*δe*(r-μy*δ)+ϵ*r*(δ-1)-1))*(b/a);
xFBe=((γ+ϵ*γe*(r-μx*γ)+ϵ*r*(γ-1))/(γ+ϵ*γe*(r-μx*γ)+ϵ*r*(γ-1)-1))*((b+β)/α);
(*Domain of PDE*)
ydagger=yFBe*daggertobarratio;
xdagger=xFBe*daggertobarratio;
Omega=Rectangle[{0,0},{ydagger,xdagger}];

(*INTERMEDIATE FUNCTIONS*)
(*Obstacles*)
Phi00[x_]:=-b + ({{\[Piecewise], ({{((((x/xstar)))^γ)*(α*xstar-β), 0<x<=xstar}, {α*x-β, x>=xstar }})}});
Phi[y_,x_]:=a*y-b + ({{\[Piecewise], ({{((((x/xstar)))^γ)*(α*xstar-β), 0<x<=xstar}, {α*x-β, x>=xstar }})}});
(*Define the function V^0(y)*)
V0[y_]:=
({{\[Piecewise], ({{b/(δ-1)*((((y/yFB)))^δ), 0<y<yFB}, {a*y-b, y>=yFB}})}});
(*Define the function V_0(x)*)
V00[x_]:=
({{\[Piecewise], ({{(α*Max[xstar,xFB]-(β+b))*((((x/(Max[xstar,xFB]))))^γ), 0<=x<Max[xstar,xFB]}, {α*x-(β+b), x>=Max[xstar,xFB]}})}});

(*DIFFERENTIAL OPERATORS*)
(*Differential operator of the ODE in y*)
L0[v_,y_]:=(1/2)((σy)^2) (y^2) D[v,{y,2}]+μy y D[v,y]-r v;
(*Differential operator of the ODE in x*)
L00[v_,x_]:=(1/2)((σx)^2) (x^2) D[v,{x,2}]+μx x D[v,x]-r v;
(*Differential operator L*)
L[v_,y_,x_]:=(1/2)((σy)^2) (y^2) D[v,{y,2}]+(1/2) ((σx)^2)(x^2) D[v,{x,2}]+ρ σy σx y x D[v,y,x]+μy y D[v,y]+μx x D[v,x]-r v;

(*BOUNDARY CONDITIONS OF THE BENCHMARK PENALIZED ODEs*)
(*Boundary conditions for the penalized ODE in y*)
bc0={
DirichletCondition[ψ[y]==b,y==0],DirichletCondition[ψ[y]==0,y==ydagger]
};
(*Boundary conditions for the penalized ODE in x*)
bc00={
DirichletCondition[ψ[x]==b,x==0],DirichletCondition[ψ[x]==0,x==xdagger]
};
(*Boundary conditions for the penalized PDE*)
bc={
DirichletCondition[ψ[y,x]==V0[y]-(a y-b),x==0],
DirichletCondition[ψ[y,x]==V00[x]-Phi00[x],y==0],
DirichletCondition[ψ[y,x]==0,x==xdagger],
DirichletCondition[ψ[y,x]==0,y==ydagger]};

(*DIFFERENTIAL EQUATIONS*)
(*Nonlinear ODE for the benchmark penalized problem in y*)
nonlinearODE0=L0[ψ[y],y]+(1/ϵ)Max[- ψ[y];0]==-L0[a*y-b,y];
(*Nonlinear ODE for the benchmark penalized problem in x*)
nonlinearODE00=L00[ψ[x],x]+(1/ϵ)Max[- ψ[x];0]==-L00[Phi00[x],x];
(*Nonlinear PDE for the two-dimensional penalized problem*)
nonlinearPDE=L[ψ[y,x],y,x]+(1/ϵ)Max[ -ψ[y,x];0]==-L[Phi[y,x],y,x];

(*NONLINEAR ODE/PDE SOLVERs*)
ψ0=NDSolveValue[
{nonlinearODE0,bc0},ψ,{y,0,ydagger},
Method->"FiniteElement"];
ψ00=NDSolveValue[
{nonlinearODE00,bc00},ψ,{x,0,xdagger},
Method->"FiniteElement"];
V00e[x_]:=ψ00[x]+Phi00[x];
ψ=NDSolveValue[
{nonlinearPDE,bc},ψ,{y,x}∈Omega,
Method->"FiniteElement"
];
];
Export["ψρ1.wxf",ψ];

Comparative statics graph for ρ

ψbaseline=Import["ψbaseline.wxf"];
xdagger=Import["xdagger.wxf"];
ydagger=Import["ydagger.wxf"];
ψ0=Import["ψρ0.wxf"];
ψ1=Import["ψρ1.wxf"];
Fig4f1=ContourPlot[ψ0[y,x]==(1/(10^4)),{x,xsmall,xdagger},{y,0,ydagger},
ContourStyle->{Directive[Black,Dotted,Thickness[0.005]]},Frame->{True,True,False,False},
FrameTicks->{{{1,2,3,4},None},{{0,50,100,150,200,250},None}},
PlotRange->{{0,xlarge},{0,ylarge}},
FrameLabel->{"Property price x","Fare revenues y","",""}];Fig4f2=ContourPlot[ψ1[y,x]==(1/(10^4)),{x,xsmall,xdagger},{y,0,ydagger},Frame->True,
ContourStyle->{Directive[Black,Dashed,Thickness[0.005]]}];
Fig4f3=ContourPlot[ψbaseline[y,x]==(1/(10^4)),{x,xsmall,xdagger},{y,0,ydagger},ContourStyle->{Directive[Black,Plain,Thickness[0.005]]}];
Figure4f=Legended[
Show[{Fig4f1,Fig4f2,Fig4f3},AspectRatio->aspectratio],Placed[Framed@LineLegend[{Directive[Black,Dotted,Thickness[0.005]],Directive[Black,Dashed,Thickness[0.005]],Directive[Black,Plain,Thickness[0.005]]},{"ρ=-.8","ρ=0","ρ=.8"},LegendMarkerSize->15,Spacings->.2,LegendLabel->"Curve (Subscript[Γ, ϵ])(·) for"],{Center,Top}]];
Export["Figure4f.png",Figure4f];SystemOpen["Figure4f.png"]

Panel G -- Negotiation power α

α=0.48 (low)

Clear[ψ,xdagger,ydagger,ψ0,ψ1];
Module[
{
(*PRIMARY PARAMETERS*)
(*Process parameters*)
σx=.2,
σy=.0675,
ρ=0.8,
μx=0.0725,
μy=0.049,
r=0.1133,
(*Profit parameters*)
I1=24.849,
c=.9,
I2=14.5,
K=4.53,
α=.48 (*CHANGED PARAMETER*)
},
(*SECONDARY PARAMETERS*)
(*Parameters for obstacle*)
a=(1/(r-μy));
b=(c/r)+I1;
β=α*I2-K;
γ=((-2 μx+((σx)^2)+2 (Sqrt[2 r ((σx)^2)+(((μx-(((σx)^2)/2)))^2)]))/(2 ((σx)^2)));
γtilde=((-2 μx+((σx)^2)-2 (Sqrt[2 r ((σx)^2)+(((μx-(((σx)^2)/2)))^2)]))/(2 ((σx)^2)));
δ=((-2 μy+((σy)^2)+2 (Sqrt[2 r ((σy)^2)+(((μy-(((σy)^2)/2)))^2)]))/(2 ((σy)^2)));
xstar=(γ/(γ-1))*(I2+(K/(1-α)));
yFB=(δ/(δ-1))*(b/a);
xFB=(γ/(γ-1))*((b+β)/α);
(*Parameters for benchmark problems*)
δe=((-2 μy+((σy)^2)+2 (Sqrt[-2 (-r-(1/ϵ)) ((σy)^2)+(((μy-(((σy)^2)/2)))^2)]))/(2 ((σy)^2)));
γe=((-2 μx+((σx)^2)+2 (Sqrt[-2 (-r-(1/ϵ)) ((σx)^2)+(((μx-(((σx)^2)/2)))^2)]))/(2 ((σx)^2)));
yFBe=((δ+ϵ*δe*(r-μy*δ)+ϵ*r*(δ-1))/(δ+ϵ*δe*(r-μy*δ)+ϵ*r*(δ-1)-1))*(b/a);
xFBe=((γ+ϵ*γe*(r-μx*γ)+ϵ*r*(γ-1))/(γ+ϵ*γe*(r-μx*γ)+ϵ*r*(γ-1)-1))*((b+β)/α);
(*Domain of PDE*)
ydagger=yFBe*daggertobarratio;
xdagger=xFBe*daggertobarratio;
Omega=Rectangle[{0,0},{ydagger,xdagger}];

(*INTERMEDIATE FUNCTIONS*)
(*Obstacles*)
Phi00[x_]:=-b + ({{\[Piecewise], ({{((((x/xstar)))^γ)*(α*xstar-β), 0<x<=xstar}, {α*x-β, x>=xstar }})}});
Phi[y_,x_]:=a*y-b + ({{\[Piecewise], ({{((((x/xstar)))^γ)*(α*xstar-β), 0<x<=xstar}, {α*x-β, x>=xstar }})}});
(*Define the function V^0(y)*)
V0[y_]:=
({{\[Piecewise], ({{b/(δ-1)*((((y/yFB)))^δ), 0<y<yFB}, {a*y-b, y>=yFB}})}});
(*Define the function V_0(x)*)
V00[x_]:=
({{\[Piecewise], ({{(α*Max[xstar,xFB]-(β+b))*((((x/(Max[xstar,xFB]))))^γ), 0<=x<Max[xstar,xFB]}, {α*x-(β+b), x>=Max[xstar,xFB]}})}});

(*DIFFERENTIAL OPERATORS*)
(*Differential operator of the ODE in y*)
L0[v_,y_]:=(1/2)((σy)^2) (y^2) D[v,{y,2}]+μy y D[v,y]-r v;
(*Differential operator of the ODE in x*)
L00[v_,x_]:=(1/2)((σx)^2) (x^2) D[v,{x,2}]+μx x D[v,x]-r v;
(*Differential operator L*)
L[v_,y_,x_]:=(1/2)((σy)^2) (y^2) D[v,{y,2}]+(1/2) ((σx)^2)(x^2) D[v,{x,2}]+ρ σy σx y x D[v,y,x]+μy y D[v,y]+μx x D[v,x]-r v;

(*BOUNDARY CONDITIONS OF THE BENCHMARK PENALIZED ODEs*)
(*Boundary conditions for the penalized ODE in y*)
bc0={
DirichletCondition[ψ[y]==b,y==0],DirichletCondition[ψ[y]==0,y==ydagger]
};
(*Boundary conditions for the penalized ODE in x*)
bc00={
DirichletCondition[ψ[x]==b,x==0],DirichletCondition[ψ[x]==0,x==xdagger]
};
(*Boundary conditions for the penalized PDE*)
bc={
DirichletCondition[ψ[y,x]==V0[y]-(a y-b),x==0],
DirichletCondition[ψ[y,x]==V00[x]-Phi00[x],y==0],
DirichletCondition[ψ[y,x]==0,x==xdagger],
DirichletCondition[ψ[y,x]==0,y==ydagger]};

(*DIFFERENTIAL EQUATIONS*)
(*Nonlinear ODE for the benchmark penalized problem in y*)
nonlinearODE0=L0[ψ[y],y]+(1/ϵ)Max[- ψ[y];0]==-L0[a*y-b,y];
(*Nonlinear ODE for the benchmark penalized problem in x*)
nonlinearODE00=L00[ψ[x],x]+(1/ϵ)Max[- ψ[x];0]==-L00[Phi00[x],x];
(*Nonlinear PDE for the two-dimensional penalized problem*)
nonlinearPDE=L[ψ[y,x],y,x]+(1/ϵ)Max[ -ψ[y,x];0]==-L[Phi[y,x],y,x];

(*NONLINEAR ODE/PDE SOLVERs*)
ψ0=NDSolveValue[
{nonlinearODE0,bc0},ψ,{y,0,ydagger},
Method->"FiniteElement"];
ψ00=NDSolveValue[
{nonlinearODE00,bc00},ψ,{x,0,xdagger},
Method->"FiniteElement"];
V00e[x_]:=ψ00[x]+Phi00[x];
ψ=NDSolveValue[
{nonlinearPDE,bc},ψ,{y,x}∈Omega,
Method->"FiniteElement"
];
];
Export["ψα0.wxf",ψ];

α=.54 (large)

Clear[ψ,xdagger,ydagger,ψ0,ψ1];
Module[
{
(*PRIMARY PARAMETERS*)
(*Process parameters*)
σx=.2,
σy=.0675,
ρ=0.8,
μx=0.0725, 
μy=0.049,
r=0.1133,
(*Profit parameters*)
I1=24.849,
c=.9,
I2=14.5,
K=4.53,
α=.54(*CHANGED PARAMETER*)
},
(*SECONDARY PARAMETERS*)
(*Parameters for obstacle*)
a=(1/(r-μy));
b=(c/r)+I1;
β=α*I2-K;
γ=((-2 μx+((σx)^2)+2 (Sqrt[2 r ((σx)^2)+(((μx-(((σx)^2)/2)))^2)]))/(2 ((σx)^2)));
γtilde=((-2 μx+((σx)^2)-2 (Sqrt[2 r ((σx)^2)+(((μx-(((σx)^2)/2)))^2)]))/(2 ((σx)^2)));
δ=((-2 μy+((σy)^2)+2 (Sqrt[2 r ((σy)^2)+(((μy-(((σy)^2)/2)))^2)]))/(2 ((σy)^2)));
xstar=(γ/(γ-1))*(I2+(K/(1-α)));
yFB=(δ/(δ-1))*(b/a);
xFB=(γ/(γ-1))*((b+β)/α);
(*Parameters for benchmark problems*)
δe=((-2 μy+((σy)^2)+2 (Sqrt[-2 (-r-(1/ϵ)) ((σy)^2)+(((μy-(((σy)^2)/2)))^2)]))/(2 ((σy)^2)));
γe=((-2 μx+((σx)^2)+2 (Sqrt[-2 (-r-(1/ϵ)) ((σx)^2)+(((μx-(((σx)^2)/2)))^2)]))/(2 ((σx)^2)));
yFBe=((δ+ϵ*δe*(r-μy*δ)+ϵ*r*(δ-1))/(δ+ϵ*δe*(r-μy*δ)+ϵ*r*(δ-1)-1))*(b/a);
xFBe=((γ+ϵ*γe*(r-μx*γ)+ϵ*r*(γ-1))/(γ+ϵ*γe*(r-μx*γ)+ϵ*r*(γ-1)-1))*((b+β)/α);
(*Domain of PDE*)
ydagger=yFBe*daggertobarratio;
xdagger=xFBe*daggertobarratio;
Omega=Rectangle[{0,0},{ydagger,xdagger}];

(*INTERMEDIATE FUNCTIONS*)
(*Obstacles*)
Phi00[x_]:=-b + ({{\[Piecewise], ({{((((x/xstar)))^γ)*(α*xstar-β), 0<x<=xstar}, {α*x-β, x>=xstar }})}});
Phi[y_,x_]:=a*y-b + ({{\[Piecewise], ({{((((x/xstar)))^γ)*(α*xstar-β), 0<x<=xstar}, {α*x-β, x>=xstar }})}});
(*Define the function V^0(y)*)
V0[y_]:=
({{\[Piecewise], ({{b/(δ-1)*((((y/yFB)))^δ), 0<y<yFB}, {a*y-b, y>=yFB}})}});
(*Define the function V_0(x)*)
V00[x_]:=
({{\[Piecewise], ({{(α*Max[xstar,xFB]-(β+b))*((((x/(Max[xstar,xFB]))))^γ), 0<=x<Max[xstar,xFB]}, {α*x-(β+b), x>=Max[xstar,xFB]}})}});

(*DIFFERENTIAL OPERATORS*)
(*Differential operator of the ODE in y*)
L0[v_,y_]:=(1/2)((σy)^2) (y^2) D[v,{y,2}]+μy y D[v,y]-r v;
(*Differential operator of the ODE in x*)
L00[v_,x_]:=(1/2)((σx)^2) (x^2) D[v,{x,2}]+μx x D[v,x]-r v;
(*Differential operator L*)
L[v_,y_,x_]:=(1/2)((σy)^2) (y^2) D[v,{y,2}]+(1/2) ((σx)^2)(x^2) D[v,{x,2}]+ρ σy σx y x D[v,y,x]+μy y D[v,y]+μx x D[v,x]-r v;

(*BOUNDARY CONDITIONS OF THE BENCHMARK PENALIZED ODEs*)
(*Boundary conditions for the penalized ODE in y*)
bc0={
DirichletCondition[ψ[y]==b,y==0],DirichletCondition[ψ[y]==0,y==ydagger]
};
(*Boundary conditions for the penalized ODE in x*)
bc00={
DirichletCondition[ψ[x]==b,x==0],DirichletCondition[ψ[x]==0,x==xdagger]
};
(*Boundary conditions for the penalized PDE*)
bc={
DirichletCondition[ψ[y,x]==V0[y]-(a y-b),x==0],
DirichletCondition[ψ[y,x]==V00[x]-Phi00[x],y==0],
DirichletCondition[ψ[y,x]==0,x==xdagger],
DirichletCondition[ψ[y,x]==0,y==ydagger]};

(*DIFFERENTIAL EQUATIONS*)
(*Nonlinear ODE for the benchmark penalized problem in y*)
nonlinearODE0=L0[ψ[y],y]+(1/ϵ)Max[- ψ[y];0]==-L0[a*y-b,y];
(*Nonlinear ODE for the benchmark penalized problem in x*)
nonlinearODE00=L00[ψ[x],x]+(1/ϵ)Max[- ψ[x];0]==-L00[Phi00[x],x];
(*Nonlinear PDE for the two-dimensional penalized problem*)
nonlinearPDE=L[ψ[y,x],y,x]+(1/ϵ)Max[ -ψ[y,x];0]==-L[Phi[y,x],y,x];

(*NONLINEAR ODE/PDE SOLVERs*)
ψ0=NDSolveValue[
{nonlinearODE0,bc0},ψ,{y,0,ydagger},
Method->"FiniteElement"];
ψ00=NDSolveValue[
{nonlinearODE00,bc00},ψ,{x,0,xdagger},
Method->"FiniteElement"];
V00e[x_]:=ψ00[x]+Phi00[x];
ψ=NDSolveValue[
{nonlinearPDE,bc},ψ,{y,x}∈Omega,
Method->"FiniteElement"
];
];
Export["ψα1.wxf",ψ];

Comparative statics graph for K

ψbaseline=Import["ψbaseline.wxf"];
xdagger=Import["xdagger.wxf"];
ydagger=Import["ydagger.wxf"];
ψ0=Import["ψα0.wxf"];
ψ1=Import["ψα1.wxf"];

Fig4g1=ContourPlot[ψ0[y,x]==(1/(10^4)),{x,xsmall,xdagger},{y,0,ydagger},
ContourStyle->{Directive[Black,Dotted,Thickness[0.005]]},
Frame->{True,True,False,False},
PlotRange->{{0,xlarge},{0,ylarge}},
FrameTicks->{{{1,2,3,4},None},{{0,50,100,150,200,250},None}},FrameLabel->{"Property price x","Fare revenues y","",""}
];
Fig4g2=ContourPlot[ψbaseline[y,x]==(1/(10^4)),{x,xsmall,xdagger},{y,0,ydagger},Frame->True,
ContourStyle->{Directive[Black,Plain,Thickness[0.005]]}];
Fig4g3=ContourPlot[ψ1[y,x]==(1/(10^4)),{x,xsmall,xdagger},{y,0,ydagger},ContourStyle->{Directive[Black,Dashed,Thickness[0.005]]}];
Figure4g=Legended[
Show[{Fig4g1,Fig4g2,Fig4g3},AspectRatio->aspectratio],Placed[Framed@LineLegend[{Directive[Black,Dotted,Thickness[0.005]],Directive[Black,Plain,Thickness[0.005]],Directive[Black,Dashed,Thickness[0.005]]},{"α=.48","α=.51","α=.54"},LegendMarkerSize->15,Spacings->.2,LegendLabel->"Curve (Subscript[Γ, ϵ])(·) for"],{Center,Top}]];
Export["Figure4g.png",Figure4g];SystemOpen["Figure4g.png"]

Panel H -- Cost K

K=0.00 (low)

Clear[ψ,xdagger,ydagger,ψ0,ψ1];
Module[
{
(*PRIMARY PARAMETERS*)
(*Process parameters*)
σx=.2,
σy=.0675,
ρ=0.8,
μx=0.0725,
μy=0.049,
r=0.1133,
(*Profit parameters*)
I1=24.849,
c=.9,
I2=14.5,
K=0, (*CHANGED PARAMETER*)
α=.51
},
(*SECONDARY PARAMETERS*)
(*Parameters for obstacle*)
a=(1/(r-μy));
b=(c/r)+I1;
β=α*I2-K;
γ=((-2 μx+((σx)^2)+2 (Sqrt[2 r ((σx)^2)+(((μx-(((σx)^2)/2)))^2)]))/(2 ((σx)^2)));
γtilde=((-2 μx+((σx)^2)-2 (Sqrt[2 r ((σx)^2)+(((μx-(((σx)^2)/2)))^2)]))/(2 ((σx)^2)));
δ=((-2 μy+((σy)^2)+2 (Sqrt[2 r ((σy)^2)+(((μy-(((σy)^2)/2)))^2)]))/(2 ((σy)^2)));
xstar=(γ/(γ-1))*(I2+(K/(1-α)));
yFB=(δ/(δ-1))*(b/a);
xFB=(γ/(γ-1))*((b+β)/α);
(*Parameters for benchmark problems*)
δe=((-2 μy+((σy)^2)+2 (Sqrt[-2 (-r-(1/ϵ)) ((σy)^2)+(((μy-(((σy)^2)/2)))^2)]))/(2 ((σy)^2)));
γe=((-2 μx+((σx)^2)+2 (Sqrt[-2 (-r-(1/ϵ)) ((σx)^2)+(((μx-(((σx)^2)/2)))^2)]))/(2 ((σx)^2)));
yFBe=((δ+ϵ*δe*(r-μy*δ)+ϵ*r*(δ-1))/(δ+ϵ*δe*(r-μy*δ)+ϵ*r*(δ-1)-1))*(b/a);
xFBe=((γ+ϵ*γe*(r-μx*γ)+ϵ*r*(γ-1))/(γ+ϵ*γe*(r-μx*γ)+ϵ*r*(γ-1)-1))*((b+β)/α);
(*Domain of PDE*)
ydagger=yFBe*daggertobarratio;
xdagger=xFBe*daggertobarratio;
Omega=Rectangle[{0,0},{ydagger,xdagger}];

(*INTERMEDIATE FUNCTIONS*)
(*Obstacles*)
Phi00[x_]:=-b + ({{\[Piecewise], ({{((((x/xstar)))^γ)*(α*xstar-β), 0<x<=xstar}, {α*x-β, x>=xstar }})}});
Phi[y_,x_]:=a*y-b + ({{\[Piecewise], ({{((((x/xstar)))^γ)*(α*xstar-β), 0<x<=xstar}, {α*x-β, x>=xstar }})}});
(*Define the function V^0(y)*)
V0[y_]:=
({{\[Piecewise], ({{b/(δ-1)*((((y/yFB)))^δ), 0<y<yFB}, {a*y-b, y>=yFB}})}});
(*Define the function V_0(x)*)
V00[x_]:=
({{\[Piecewise], ({{(α*Max[xstar,xFB]-(β+b))*((((x/(Max[xstar,xFB]))))^γ), 0<=x<Max[xstar,xFB]}, {α*x-(β+b), x>=Max[xstar,xFB]}})}});

(*DIFFERENTIAL OPERATORS*)
(*Differential operator of the ODE in y*)
L0[v_,y_]:=(1/2)((σy)^2) (y^2) D[v,{y,2}]+μy y D[v,y]-r v;
(*Differential operator of the ODE in x*)
L00[v_,x_]:=(1/2)((σx)^2) (x^2) D[v,{x,2}]+μx x D[v,x]-r v;
(*Differential operator L*)
L[v_,y_,x_]:=(1/2)((σy)^2) (y^2) D[v,{y,2}]+(1/2) ((σx)^2)(x^2) D[v,{x,2}]+ρ σy σx y x D[v,y,x]+μy y D[v,y]+μx x D[v,x]-r v;

(*BOUNDARY CONDITIONS OF THE BENCHMARK PENALIZED ODEs*)
(*Boundary conditions for the penalized ODE in y*)
bc0={
DirichletCondition[ψ[y]==b,y==0],DirichletCondition[ψ[y]==0,y==ydagger]
};
(*Boundary conditions for the penalized ODE in x*)
bc00={
DirichletCondition[ψ[x]==b,x==0],DirichletCondition[ψ[x]==0,x==xdagger]
};
(*Boundary conditions for the penalized PDE*)
bc={
DirichletCondition[ψ[y,x]==V0[y]-(a y-b),x==0],
DirichletCondition[ψ[y,x]==V00[x]-Phi00[x],y==0],
DirichletCondition[ψ[y,x]==0,x==xdagger],
DirichletCondition[ψ[y,x]==0,y==ydagger]};

(*DIFFERENTIAL EQUATIONS*)
(*Nonlinear ODE for the benchmark penalized problem in y*)
nonlinearODE0=L0[ψ[y],y]+(1/ϵ)Max[- ψ[y];0]==-L0[a*y-b,y];
(*Nonlinear ODE for the benchmark penalized problem in x*)
nonlinearODE00=L00[ψ[x],x]+(1/ϵ)Max[- ψ[x];0]==-L00[Phi00[x],x];
(*Nonlinear PDE for the two-dimensional penalized problem*)
nonlinearPDE=L[ψ[y,x],y,x]+(1/ϵ)Max[ -ψ[y,x];0]==-L[Phi[y,x],y,x];

(*NONLINEAR ODE/PDE SOLVERs*)
ψ0=NDSolveValue[
{nonlinearODE0,bc0},ψ,{y,0,ydagger},
Method->"FiniteElement"];
ψ00=NDSolveValue[
{nonlinearODE00,bc00},ψ,{x,0,xdagger},
Method->"FiniteElement"];
V00e[x_]:=ψ00[x]+Phi00[x];
ψ=NDSolveValue[
{nonlinearPDE,bc},ψ,{y,x}∈Omega,
Method->"FiniteElement"
];
];
Export["ψK0.wxf",ψ];

K=9.00 (large)

Clear[ψ,xdagger,ydagger,ψ0,ψ1];
Module[
{
(*PRIMARY PARAMETERS*)
(*Process parameters*)
σx=.2,
σy=.0675,
ρ=0.8,
μx=0.0725, 
μy=0.049,
r=0.1133,
(*Profit parameters*)
I1=24.849,
c=.9,
I2=14.5,
K=9.00,(*CHANGED PARAMETER*)
α=.51
},
(*SECONDARY PARAMETERS*)
(*Parameters for obstacle*)
a=(1/(r-μy));
b=(c/r)+I1;
β=α*I2-K;
γ=((-2 μx+((σx)^2)+2 (Sqrt[2 r ((σx)^2)+(((μx-(((σx)^2)/2)))^2)]))/(2 ((σx)^2)));
γtilde=((-2 μx+((σx)^2)-2 (Sqrt[2 r ((σx)^2)+(((μx-(((σx)^2)/2)))^2)]))/(2 ((σx)^2)));
δ=((-2 μy+((σy)^2)+2 (Sqrt[2 r ((σy)^2)+(((μy-(((σy)^2)/2)))^2)]))/(2 ((σy)^2)));
xstar=(γ/(γ-1))*(I2+(K/(1-α)));
yFB=(δ/(δ-1))*(b/a);
xFB=(γ/(γ-1))*((b+β)/α);
(*Parameters for benchmark problems*)
δe=((-2 μy+((σy)^2)+2 (Sqrt[-2 (-r-(1/ϵ)) ((σy)^2)+(((μy-(((σy)^2)/2)))^2)]))/(2 ((σy)^2)));
γe=((-2 μx+((σx)^2)+2 (Sqrt[-2 (-r-(1/ϵ)) ((σx)^2)+(((μx-(((σx)^2)/2)))^2)]))/(2 ((σx)^2)));
yFBe=((δ+ϵ*δe*(r-μy*δ)+ϵ*r*(δ-1))/(δ+ϵ*δe*(r-μy*δ)+ϵ*r*(δ-1)-1))*(b/a);
xFBe=((γ+ϵ*γe*(r-μx*γ)+ϵ*r*(γ-1))/(γ+ϵ*γe*(r-μx*γ)+ϵ*r*(γ-1)-1))*((b+β)/α);
(*Domain of PDE*)
ydagger=yFBe*daggertobarratio;
xdagger=xFBe*daggertobarratio;
Omega=Rectangle[{0,0},{ydagger,xdagger}];

(*INTERMEDIATE FUNCTIONS*)
(*Obstacles*)
Phi00[x_]:=-b + ({{\[Piecewise], ({{((((x/xstar)))^γ)*(α*xstar-β), 0<x<=xstar}, {α*x-β, x>=xstar }})}});
Phi[y_,x_]:=a*y-b + ({{\[Piecewise], ({{((((x/xstar)))^γ)*(α*xstar-β), 0<x<=xstar}, {α*x-β, x>=xstar }})}});
(*Define the function V^0(y)*)
V0[y_]:=
({{\[Piecewise], ({{b/(δ-1)*((((y/yFB)))^δ), 0<y<yFB}, {a*y-b, y>=yFB}})}});
(*Define the function V_0(x)*)
V00[x_]:=
({{\[Piecewise], ({{(α*Max[xstar,xFB]-(β+b))*((((x/(Max[xstar,xFB]))))^γ), 0<=x<Max[xstar,xFB]}, {α*x-(β+b), x>=Max[xstar,xFB]}})}});

(*DIFFERENTIAL OPERATORS*)
(*Differential operator of the ODE in y*)
L0[v_,y_]:=(1/2)((σy)^2) (y^2) D[v,{y,2}]+μy y D[v,y]-r v;
(*Differential operator of the ODE in x*)
L00[v_,x_]:=(1/2)((σx)^2) (x^2) D[v,{x,2}]+μx x D[v,x]-r v;
(*Differential operator L*)
L[v_,y_,x_]:=(1/2)((σy)^2) (y^2) D[v,{y,2}]+(1/2) ((σx)^2)(x^2) D[v,{x,2}]+ρ σy σx y x D[v,y,x]+μy y D[v,y]+μx x D[v,x]-r v;

(*BOUNDARY CONDITIONS OF THE BENCHMARK PENALIZED ODEs*)
(*Boundary conditions for the penalized ODE in y*)
bc0={
DirichletCondition[ψ[y]==b,y==0],DirichletCondition[ψ[y]==0,y==ydagger]
};
(*Boundary conditions for the penalized ODE in x*)
bc00={
DirichletCondition[ψ[x]==b,x==0],DirichletCondition[ψ[x]==0,x==xdagger]
};
(*Boundary conditions for the penalized PDE*)
bc={
DirichletCondition[ψ[y,x]==V0[y]-(a y-b),x==0],
DirichletCondition[ψ[y,x]==V00[x]-Phi00[x],y==0],
DirichletCondition[ψ[y,x]==0,x==xdagger],
DirichletCondition[ψ[y,x]==0,y==ydagger]};

(*DIFFERENTIAL EQUATIONS*)
(*Nonlinear ODE for the benchmark penalized problem in y*)
nonlinearODE0=L0[ψ[y],y]+(1/ϵ)Max[- ψ[y];0]==-L0[a*y-b,y];
(*Nonlinear ODE for the benchmark penalized problem in x*)
nonlinearODE00=L00[ψ[x],x]+(1/ϵ)Max[- ψ[x];0]==-L00[Phi00[x],x];
(*Nonlinear PDE for the two-dimensional penalized problem*)
nonlinearPDE=L[ψ[y,x],y,x]+(1/ϵ)Max[ -ψ[y,x];0]==-L[Phi[y,x],y,x];

(*NONLINEAR ODE/PDE SOLVERs*)
ψ0=NDSolveValue[
{nonlinearODE0,bc0},ψ,{y,0,ydagger},
Method->"FiniteElement"];
ψ00=NDSolveValue[
{nonlinearODE00,bc00},ψ,{x,0,xdagger},
Method->"FiniteElement"];
V00e[x_]:=ψ00[x]+Phi00[x];
ψ=NDSolveValue[
{nonlinearPDE,bc},ψ,{y,x}∈Omega,
Method->"FiniteElement"
];
];
Export["ψK1.wxf",ψ];

Comparative statics graph for K

ψbaseline=Import["ψbaseline.wxf"];
xdagger=Import["xdagger.wxf"];
ydagger=Import["ydagger.wxf"];
ψ0=Import["ψK0.wxf"];
ψ1=Import["ψK1.wxf"];

Fig4h1=ContourPlot[ψ0[y,x]==(1/(10^4)),{x,xsmall,xdagger},{y,0,ydagger},
ContourStyle->{Directive[Black,Dotted,Thickness[0.005]]},
Frame->{True,True,False,False},
PlotRange->{{0,xlarge},{0,ylarge}},
FrameTicks->{{{1,2,3,4},None},{{0,50,100,150,200,250},None}},FrameLabel->{"Property price x","Fare revenues y","",""}
];
Fig4h2=ContourPlot[ψbaseline[y,x]==(1/(10^4)),{x,xsmall,xdagger},{y,0,ydagger},Frame->True,
ContourStyle->{Directive[Black,Plain,Thickness[0.005]]}];
Fig4h3=ContourPlot[ψ1[y,x]==(1/(10^4)),{x,xsmall,xdagger},{y,0,ydagger},ContourStyle->{Directive[Black,Dashed,Thickness[0.005]]}];
Figure4h=Legended[
Show[{Fig4h1,Fig4h2,Fig4h3},AspectRatio->aspectratio],Placed[Framed@LineLegend[{Directive[Black,Dotted,Thickness[0.005]],Directive[Black,Dashed,Thickness[0.005]],Directive[Black,Plain,Thickness[0.005]]},{"K=.00","K=4.53","K=9.00"},LegendMarkerSize->15,Spacings->.2,LegendLabel->"Curve (Subscript[Γ, ϵ])(·) for"],{Center,Top}]];
Export["Figure4h.png",Figure4h];SystemOpen["Figure4h.png"]

Panel I -- Cost (Subscript[I, 2])

(Subscript[I, 2])=0.0 (low)

Clear[ψ,xdagger,ydagger,ψ0,ψ1];
Module[
{
(*PRIMARY PARAMETERS*)
(*Process parameters*)
σx=.2,
σy=.0675,
ρ=0.8,
μx=0.0725,
μy=0.049,
r=0.1133,
(*Profit parameters*)
I1=24.849,
c=.9,
I2=0,(*CHANGED PARAMETER*)
K=4.53, 
α=.51
},
(*SECONDARY PARAMETERS*)
(*Parameters for obstacle*)
a=(1/(r-μy));
b=(c/r)+I1;
β=α*I2-K;
γ=((-2 μx+((σx)^2)+2 (Sqrt[2 r ((σx)^2)+(((μx-(((σx)^2)/2)))^2)]))/(2 ((σx)^2)));
γtilde=((-2 μx+((σx)^2)-2 (Sqrt[2 r ((σx)^2)+(((μx-(((σx)^2)/2)))^2)]))/(2 ((σx)^2)));
δ=((-2 μy+((σy)^2)+2 (Sqrt[2 r ((σy)^2)+(((μy-(((σy)^2)/2)))^2)]))/(2 ((σy)^2)));
xstar=(γ/(γ-1))*(I2+(K/(1-α)));
yFB=(δ/(δ-1))*(b/a);
xFB=(γ/(γ-1))*((b+β)/α);
(*Parameters for benchmark problems*)
δe=((-2 μy+((σy)^2)+2 (Sqrt[-2 (-r-(1/ϵ)) ((σy)^2)+(((μy-(((σy)^2)/2)))^2)]))/(2 ((σy)^2)));
γe=((-2 μx+((σx)^2)+2 (Sqrt[-2 (-r-(1/ϵ)) ((σx)^2)+(((μx-(((σx)^2)/2)))^2)]))/(2 ((σx)^2)));
yFBe=((δ+ϵ*δe*(r-μy*δ)+ϵ*r*(δ-1))/(δ+ϵ*δe*(r-μy*δ)+ϵ*r*(δ-1)-1))*(b/a);
xFBe=((γ+ϵ*γe*(r-μx*γ)+ϵ*r*(γ-1))/(γ+ϵ*γe*(r-μx*γ)+ϵ*r*(γ-1)-1))*((b+β)/α);
(*Domain of PDE*)
ydagger=yFBe*daggertobarratio;
xdagger=xFBe*daggertobarratio;
Omega=Rectangle[{0,0},{ydagger,xdagger}];

(*INTERMEDIATE FUNCTIONS*)
(*Obstacles*)
Phi00[x_]:=-b + ({{\[Piecewise], ({{((((x/xstar)))^γ)*(α*xstar-β), 0<x<=xstar}, {α*x-β, x>=xstar }})}});
Phi[y_,x_]:=a*y-b + ({{\[Piecewise], ({{((((x/xstar)))^γ)*(α*xstar-β), 0<x<=xstar}, {α*x-β, x>=xstar }})}});
(*Define the function V^0(y)*)
V0[y_]:=
({{\[Piecewise], ({{b/(δ-1)*((((y/yFB)))^δ), 0<y<yFB}, {a*y-b, y>=yFB}})}});
(*Define the function V_0(x)*)
V00[x_]:=
({{\[Piecewise], ({{(α*Max[xstar,xFB]-(β+b))*((((x/(Max[xstar,xFB]))))^γ), 0<=x<Max[xstar,xFB]}, {α*x-(β+b), x>=Max[xstar,xFB]}})}});

(*DIFFERENTIAL OPERATORS*)
(*Differential operator of the ODE in y*)
L0[v_,y_]:=(1/2)((σy)^2) (y^2) D[v,{y,2}]+μy y D[v,y]-r v;
(*Differential operator of the ODE in x*)
L00[v_,x_]:=(1/2)((σx)^2) (x^2) D[v,{x,2}]+μx x D[v,x]-r v;
(*Differential operator L*)
L[v_,y_,x_]:=(1/2)((σy)^2) (y^2) D[v,{y,2}]+(1/2) ((σx)^2)(x^2) D[v,{x,2}]+ρ σy σx y x D[v,y,x]+μy y D[v,y]+μx x D[v,x]-r v;

(*BOUNDARY CONDITIONS OF THE BENCHMARK PENALIZED ODEs*)
(*Boundary conditions for the penalized ODE in y*)
bc0={
DirichletCondition[ψ[y]==b,y==0],DirichletCondition[ψ[y]==0,y==ydagger]
};
(*Boundary conditions for the penalized ODE in x*)
bc00={
DirichletCondition[ψ[x]==b,x==0],DirichletCondition[ψ[x]==0,x==xdagger]
};
(*Boundary conditions for the penalized PDE*)
bc={
DirichletCondition[ψ[y,x]==V0[y]-(a y-b),x==0],
DirichletCondition[ψ[y,x]==V00[x]-Phi00[x],y==0],
DirichletCondition[ψ[y,x]==0,x==xdagger],
DirichletCondition[ψ[y,x]==0,y==ydagger]};

(*DIFFERENTIAL EQUATIONS*)
(*Nonlinear ODE for the benchmark penalized problem in y*)
nonlinearODE0=L0[ψ[y],y]+(1/ϵ)Max[- ψ[y];0]==-L0[a*y-b,y];
(*Nonlinear ODE for the benchmark penalized problem in x*)
nonlinearODE00=L00[ψ[x],x]+(1/ϵ)Max[- ψ[x];0]==-L00[Phi00[x],x];
(*Nonlinear PDE for the two-dimensional penalized problem*)
nonlinearPDE=L[ψ[y,x],y,x]+(1/ϵ)Max[ -ψ[y,x];0]==-L[Phi[y,x],y,x];

(*NONLINEAR ODE/PDE SOLVERs*)
ψ0=NDSolveValue[
{nonlinearODE0,bc0},ψ,{y,0,ydagger},
Method->"FiniteElement"];
ψ00=NDSolveValue[
{nonlinearODE00,bc00},ψ,{x,0,xdagger},
Method->"FiniteElement"];
V00e[x_]:=ψ00[x]+Phi00[x];
ψ=NDSolveValue[
{nonlinearPDE,bc},ψ,{y,x}∈Omega,
Method->"FiniteElement"
];
];
Export["ψI0.wxf",ψ];

(Subscript[I, 2])=29.0 (large)

Clear[ψ,xdagger,ydagger,ψ0,ψ1];
Module[
{
(*PRIMARY PARAMETERS*)
(*Process parameters*)
σx=.2,
σy=.0675,
ρ=0.8,
μx=0.0725, 
μy=0.049,
r=0.1133,
(*Profit parameters*)
I1=24.849,
c=.9,
I2=29,(*CHANGED PARAMETER*)
K=4.53,
α=.51
},
(*SECONDARY PARAMETERS*)
(*Parameters for obstacle*)
a=(1/(r-μy));
b=(c/r)+I1;
β=α*I2-K;
γ=((-2 μx+((σx)^2)+2 (Sqrt[2 r ((σx)^2)+(((μx-(((σx)^2)/2)))^2)]))/(2 ((σx)^2)));
γtilde=((-2 μx+((σx)^2)-2 (Sqrt[2 r ((σx)^2)+(((μx-(((σx)^2)/2)))^2)]))/(2 ((σx)^2)));
δ=((-2 μy+((σy)^2)+2 (Sqrt[2 r ((σy)^2)+(((μy-(((σy)^2)/2)))^2)]))/(2 ((σy)^2)));
xstar=(γ/(γ-1))*(I2+(K/(1-α)));
yFB=(δ/(δ-1))*(b/a);
xFB=(γ/(γ-1))*((b+β)/α);
(*Parameters for benchmark problems*)
δe=((-2 μy+((σy)^2)+2 (Sqrt[-2 (-r-(1/ϵ)) ((σy)^2)+(((μy-(((σy)^2)/2)))^2)]))/(2 ((σy)^2)));
γe=((-2 μx+((σx)^2)+2 (Sqrt[-2 (-r-(1/ϵ)) ((σx)^2)+(((μx-(((σx)^2)/2)))^2)]))/(2 ((σx)^2)));
yFBe=((δ+ϵ*δe*(r-μy*δ)+ϵ*r*(δ-1))/(δ+ϵ*δe*(r-μy*δ)+ϵ*r*(δ-1)-1))*(b/a);
xFBe=((γ+ϵ*γe*(r-μx*γ)+ϵ*r*(γ-1))/(γ+ϵ*γe*(r-μx*γ)+ϵ*r*(γ-1)-1))*((b+β)/α);
(*Domain of PDE*)
ydagger=yFBe*daggertobarratio;
xdagger=xFBe*daggertobarratio;
Omega=Rectangle[{0,0},{ydagger,xdagger}];

(*INTERMEDIATE FUNCTIONS*)
(*Obstacles*)
Phi00[x_]:=-b + ({{\[Piecewise], ({{((((x/xstar)))^γ)*(α*xstar-β), 0<x<=xstar}, {α*x-β, x>=xstar }})}});
Phi[y_,x_]:=a*y-b + ({{\[Piecewise], ({{((((x/xstar)))^γ)*(α*xstar-β), 0<x<=xstar}, {α*x-β, x>=xstar }})}});
(*Define the function V^0(y)*)
V0[y_]:=
({{\[Piecewise], ({{b/(δ-1)*((((y/yFB)))^δ), 0<y<yFB}, {a*y-b, y>=yFB}})}});
(*Define the function V_0(x)*)
V00[x_]:=
({{\[Piecewise], ({{(α*Max[xstar,xFB]-(β+b))*((((x/(Max[xstar,xFB]))))^γ), 0<=x<Max[xstar,xFB]}, {α*x-(β+b), x>=Max[xstar,xFB]}})}});

(*DIFFERENTIAL OPERATORS*)
(*Differential operator of the ODE in y*)
L0[v_,y_]:=(1/2)((σy)^2) (y^2) D[v,{y,2}]+μy y D[v,y]-r v;
(*Differential operator of the ODE in x*)
L00[v_,x_]:=(1/2)((σx)^2) (x^2) D[v,{x,2}]+μx x D[v,x]-r v;
(*Differential operator L*)
L[v_,y_,x_]:=(1/2)((σy)^2) (y^2) D[v,{y,2}]+(1/2) ((σx)^2)(x^2) D[v,{x,2}]+ρ σy σx y x D[v,y,x]+μy y D[v,y]+μx x D[v,x]-r v;

(*BOUNDARY CONDITIONS OF THE BENCHMARK PENALIZED ODEs*)
(*Boundary conditions for the penalized ODE in y*)
bc0={
DirichletCondition[ψ[y]==b,y==0],DirichletCondition[ψ[y]==0,y==ydagger]
};
(*Boundary conditions for the penalized ODE in x*)
bc00={
DirichletCondition[ψ[x]==b,x==0],DirichletCondition[ψ[x]==0,x==xdagger]
};
(*Boundary conditions for the penalized PDE*)
bc={
DirichletCondition[ψ[y,x]==V0[y]-(a y-b),x==0],
DirichletCondition[ψ[y,x]==V00[x]-Phi00[x],y==0],
DirichletCondition[ψ[y,x]==0,x==xdagger],
DirichletCondition[ψ[y,x]==0,y==ydagger]};

(*DIFFERENTIAL EQUATIONS*)
(*Nonlinear ODE for the benchmark penalized problem in y*)
nonlinearODE0=L0[ψ[y],y]+(1/ϵ)Max[- ψ[y];0]==-L0[a*y-b,y];
(*Nonlinear ODE for the benchmark penalized problem in x*)
nonlinearODE00=L00[ψ[x],x]+(1/ϵ)Max[- ψ[x];0]==-L00[Phi00[x],x];
(*Nonlinear PDE for the two-dimensional penalized problem*)
nonlinearPDE=L[ψ[y,x],y,x]+(1/ϵ)Max[ -ψ[y,x];0]==-L[Phi[y,x],y,x];

(*NONLINEAR ODE/PDE SOLVERs*)
ψ0=NDSolveValue[
{nonlinearODE0,bc0},ψ,{y,0,ydagger},
Method->"FiniteElement"];
ψ00=NDSolveValue[
{nonlinearODE00,bc00},ψ,{x,0,xdagger},
Method->"FiniteElement"];
V00e[x_]:=ψ00[x]+Phi00[x];
ψ=NDSolveValue[
{nonlinearPDE,bc},ψ,{y,x}∈Omega,
Method->"FiniteElement"
];
];
Export["ψI1.wxf",ψ];

Comparative statics graph for K

ψbaseline=Import["ψbaseline.wxf"];
xdagger=Import["xdagger.wxf"];
ydagger=Import["ydagger.wxf"];
ψ0=Import["ψI0.wxf"];
ψ1=Import["ψI1.wxf"];

Fig4i1=ContourPlot[ψ0[y,x]==(1/(10^4)),{x,xsmall,xdagger},{y,0,ydagger},
ContourStyle->{Directive[Black,Dotted,Thickness[0.005]]},
Frame->{True,True,False,False},
PlotRange->{{0,xlarge},{0,ylarge}},
FrameTicks->{{{1,2,3,4},None},{{0,50,100,150,200,250},None}},FrameLabel->{"Property price x","Fare revenues y","",""}
];
Fig4i2=ContourPlot[ψbaseline[y,x]==(1/(10^4)),{x,xsmall,xdagger},{y,0,ydagger},Frame->True,
ContourStyle->{Directive[Black,Plain,Thickness[0.005]]}];
Fig4i3=ContourPlot[ψ1[y,x]==(1/(10^4)),{x,xsmall,xdagger},{y,0,ydagger},ContourStyle->{Directive[Black,Dashed,Thickness[0.005]]}];
Figure4i=Legended[
Show[{Fig4i1,Fig4i2,Fig4i3},AspectRatio->aspectratio],Placed[Framed@LineLegend[{Directive[Black,Dotted,Thickness[0.005]],Directive[Black,Plain,Thickness[0.005]],Directive[Black,Dashed,Thickness[0.005]]},{"(Subscript[I, 2])=0.0","(Subscript[I, 2])=14.5","(Subscript[I, 2])=29.0"},LegendMarkerSize->15,Spacings->.2,LegendLabel->"Curve (Subscript[Γ, ϵ])(·) for"],{Center,Top}]];
Export["Figure4i.png",Figure4i];SystemOpen["Figure4i.png"]