The net worth trap: Investment and output dynamics in the presence of financing constraints

J. Isohätälä, A. Milne and D. Robertson

Introduction / Usage

First Group of Cells with documentation and setting up parameters and solution functions.

The central functions in this notebook are SolveandEverything[] and SolveandEverythingNoRent[] . 

In order to solve for recapitalisation they allow for the maximum level of borrowing to be different from the value derived in Appendix A of the paper. The latter approach creates programming complications, because then the inititial point is not at the “singular intersection” (in the case of the model with rental this is where the analytic power law approximations for q and psi are applicable) and there are two solutions for etastar and q(etastar) that satisfy the boundary conditions, a “lower branch” and an “upper branch”. We need to avoid the upper branch because this is associated with negative dividends.   

Within the “solveeverything” functions,two “numerical solvers” are used, AuxSolveqWfF and AuxSolveqWfFnorent, which take parameter values and initial value of qI on lower bound, and spit out q, W, f and F as peice wise numerical function ( first computing q from lowest value etabar upto etastar , then computing W, f, F . ) , together with four values of eta ηI (lowest value),  ηC (where the “cut” occurs moving from power law to numerical solution, ηT (transition from renting to no rental), and ηS  the upper boundary *)

The choice between liquidation and recapitalisation is determined within the “solveeverything” functions by search carried out by functions “ComputeLowerBoundaryq” and “ComputeLowerBoundaryqNoRent”, these search on qI to enforce matching condition with value gain or liquidation. 

A further function “ComputeScaledMedian” and “ComputeScaledMedianNoRent” are used to prepare the contour charts.

Most output is graphical. For the main figures this is controlled using further functions SolveandPlotEverything[] and SoveandPlotEverythingNoRent[]. The countour charts are produced directly using Mathematica graphics commands.

The notebook contains four groups of cells. First group contains this documentation and sets up all parameters and functions. This must be run intially before doing anything else. Second “group” (a single cell provides the key graphic for rental, allowing the user to control output using sliders. )  Third group allows substittuion of parameters and then creates the Figures from the paper. 

For output purposes, preferable to run with messages window open.

Short usage:
1. See input parameters and output configuration in “Parameters and configuration” section
2. Choose from menu Evaluation -> Evaluate Notebook
3. Wait about 10 mins
4. All figures for the paper are generated, and you are done.

Or:
1. See input parameters and output configuration in “Parameters and configuration” section
2. Choose from menu Evaluation -> Evaluate Initialization Cells. This only sets up the solver algorithms
3. Run individually the sections for each plot, or play with the interactive part.

The actual solver routines and supporting functions are under their respective sections. These are “folded” away from sight in the current save of this notebook; just double click the bracket on the right-hand side to open the sections to see the internals.

Initialization

This section is marked as initialization cells. To run all cells here, just choose Evaluation -> Evaluate Initialization Cells from the menu.

Parameters and configuration

This first set of cells sets baseline parameters, from the paper, used in the numerical solution. The output path for the figures is specified here.

All graphics will be exported to the following folder PathToGraphics, extra information (typically not needed) is printed PrintExtraInfo is set to True.

PathToGraphics="C:\\Users\\Alistair\\Documents\\research\\qmodel\\graphics\\";
PathToGraphics="/home/jisohata/Dropbox/Projects/RALboro/Final/";
PathToGraphics="C:\\Users\\Jukka\ Isohatala\\Desktop\\";
PrintExtraInfo=False;

The following cell sets the default parameters, these baseline parameterisations are as follows:

({{ρ, insider discount rate, 0.06}, {r, outsider discount rate, 0.05}, {σ1, volatility of cash flow, 0.2}, {σ2, volatility of capital productivity, 0.0}, {a, productivity of insider capital, 0.1}, {ab, productivity of outsider capital, 0.04}, {θ, costs of capital adjustment, 15.0}, {δ, depreciation rate, 0.02}, {χ, costs of recapitalisation, 0.75}, {μ, Cost of adjustment parameter, 0}, {ϵCut, numerical solution parameter, 1.0 (10^(-3))}, {qBarSearchNDiv, numerical solution parameter, 10}})

These are stored in two lists “defaultParametersRule” and (for comparison purposes) “defPars”

defaultParametersRule={
ρ->0.06,r->0.05,
σ1->0.2,σ2->0.0,
a->0.1,ab->0.04,
θ->15.0,δ->0.02,
μ->0.0,
χ->0.75,
ϵCut->1.0 (10^(-3)),
qBarSearchNDiv->10
};
defPars=defaultParametersRule;

This ends the configuration / parameter initialization.
Some tunable parameters are found in the cells generating the plots, these mostly control details of the output graphics.

Auxilliary routines

Generic

(* we also have a function "overrideDefaultParameters[]" which substitutes parameter values in defaultParametersRule , example of usage ComputeScaledMedian[overrideDefaultParameters[{σ1->0.35,χ->100}]] *)
overrideDefaultParameters[rules_]:=Module[{v=(rules/.((A_->B_):>A))},
Select[defaultParametersRule,(!MemberQ[v,#1/.((A_->B_):>A)]&)]~Join~rules
];
(* This following function is simply used to dump the default values of some quantities of interest onto the screen; now disabled *)
PrintDefValue[symbol_]:=If[PrintExtraInfo,Print[symbol," = ",Symbol@symbol/.defaultParametersRule]];

Plotting and graphics

Contour and 3D plots tend to generate massive pdf files. For that reason, here are variant versions of ContourPlot, Plot3D, and ParametricPlot3D. They automatically produce rasterized output which will be far smaller in size than the equivalent vector graphics version. Here, only the plotted surface / contour is rasterized. The axes, labels, etc. are still vector graphics: This compromise should give the good quality and small size.

Clear[RasterContourPlot]
SetAttributes[RasterContourPlot,HoldAll]
RasterContourPlot[expr_,{x_,x0_,x1_},{y_,y0_,y1_},options:OptionsPattern[]]:=Module[{rplot},
rplot=ContourPlot[expr,{x,x0,x1},{y,y0,y1}
,Frame->False
,Axes->False
,Evaluate@FilterRules[{options},Options[ContourPlot]]
];
ContourPlot[I,{x,x0,x1},{y,y0,y1}
,Prolog->Inset[Rasterize[rplot,"Graphics",FilterRules[{options},Options[Rasterize]]],Scaled@{0,0},Scaled@{0,0},Scaled@{1,1}]
,Evaluate@FilterRules[{options},Options[ContourPlot]]
]
];
Clear[RasterPlot3D]
SetAttributes[RasterPlot3D,HoldAll]
RasterPlot3D[expr_,{x_,x0_,x1_},{y_,y0_,y1_},options:OptionsPattern[]]:=Module[{rplot},
rplot=Plot3D[expr,{x,x0,x1},{y,y0,y1}
,Axes->False
,Boxed->True
,Evaluate@FilterRules[{options},Options[Plot3D]]
];
Plot3D[I,{x,x0,x1},{y,y0,y1}
,Prolog->Inset[Rasterize[rplot,"Graphics",FilterRules[{options},Options[Rasterize]]],Scaled@{0,0},Scaled@{0,0},Scaled@{1,1}]
,Evaluate@FilterRules[{options},Options[Plot3D]]
]
];
Clear[RasterParametricPlot3D]
SetAttributes[RasterParametricPlot3D,HoldAll]
RasterParametricPlot3D[expr_,{x_,x0_,x1_},{y_,y0_,y1_},options:OptionsPattern[]]:=Module[{rplot},
rplot=ParametricPlot3D[expr,{x,x0,x1},{y,y0,y1}
,Axes->False
,Boxed->True
,Evaluate@FilterRules[{options},Options[Plot3D]]
];
ParametricPlot3D[I,{x,x0,x1},{y,y0,y1}
,Prolog->Inset[Rasterize[rplot,"Graphics",FilterRules[{options},Options[Rasterize]]],Scaled@{0,0},Scaled@{0,0},Scaled@{1,1}]
,Evaluate@FilterRules[{options},Options[Plot3D]]
]
];

Solver algorithm

Model with rental

(* We now introduce some functions and expressions from the paper. 
(1) In the model with rent, Ψ , the proportion of retained capital, is function of η and q as in paper (Nov 2012) eqn (24) bottom of page 28. These are the same!! though not immediately clear (possibly the version here has been derived using Mathematica  *)
Ψ[η_,q_]:=-((1+(q^2)+2 θ (ab-δ+r η+μ-η ρ)-2 q (1+θ (μ+ρ)))/((a-ab) θ));
(* (2) An equation for lower value of η; this is equation (11) page 18, of Nov 2012 version, the point at which the valuation of capital by outsiders equals the level of debt. in more general models the minimum value ηbar could be lower or higher than this (i.e. this coudl be an additional parameter) but we restrict ourselves to the case ηbar=ηsingular *)
ηSingular=-1-θ (r+μ)+(Sqrt[θ (-2 ab+2 δ+r (2+r θ)+2 r θ μ+θ (μ^2))]);
PrintDefValue["ηSingular"]
(* This should evaluate to a negative value using default parameters, it is negative because it is a debt (+ve would be the holding of a cash balance)*)

(* (3) ODE in q for case with rental. Corresponds to equation (23), page 28: 
suggestion will be easier to read if written using Ψ[η_,q_] so becomes qPrimeRentingEq={q'[η]==((a-ab)  (η+q[η]))/((((σ1)^2)+(η^2) ((σ2)^2)) Ψ[η,q])};
 *)
qPrimeRentingEq={q'[η]==-((((a-ab))^2) θ (η+q[η]))/((((σ1)^2)+(η^2) ((σ2)^2)) (1+2 θ (ab-δ+r η+μ-η ρ)-2 (1+θ (μ+ρ)) q[η]+((q[η])^2)))};
(* (4) various function for q' (numerical solution of q' is what is required to pin down the model), taking account of possible singularity if η+q=0 and ψ=0;when renting out of capital stock takes place (ψ<1) then this is equation (23) on page 28 of the Nov 2012 version of the paper; when there is no renting out of capital stock then this is equation (8) on page 17 of the Nov 2012 equation for 
 alternative expression used when there is a singularity is equation (34) from Appendix B, page 46, Nov 2012 paper 
allowing for rental and for non-rental *)
qPrimeRentingFun[η_,q_]=With[{ψ=Ψ[η,q],c=(a-ab)/(((σ1)^2)+(η^2)((σ2)^2)),γ=(1/2)(((a-ab))^2)/(((σ1)^2)+(η^2)((σ2)^2))},
If[η+q==0&&Abs[ψ]<√$MachineEpsilon,
((-(ρ-r-γ)+(Sqrt[(((ρ-r-γ))^2)+((4γ)/θ)(1+θ(ρ+μ)-q)]))/(2(θ^(-1))(1+θ(ρ+μ)-q))),
(1/ψ)c(q+η)
]
];
qPrimeNotRentingFun[η_,q_]=(((q+η) (1+(q^2)+2 θ (a-δ+r η+μ-η ρ)-2 q (1+θ (μ+ρ))))/(θ (((σ1)^2)+(η^2) ((σ2)^2))));

qPrimeRentingNearSingularEq={q'[η]==qPrimeRentingFun[η,q[η]]};
qPrimeNotRentingEq={q'[η]==qPrimeNotRentingFun[η,q[η]]};

(* Expressions for volatility, in two cases (renting, not renting)*)
halfσηSquare[η_,q_]=(1/2)(((σ1)^2)+(η^2)((σ2)^2))((Ψ[η,q])^2);
halfσηSquareNotRenting[η_,q_]=(1/2)(((σ1)^2)+(η^2)((σ2)^2));
(*The following are the equations for ergodic density, equations (26a) and (26b) on page 30 of the Nov 2012 paper (with additional parameter arising because of volatility in cap productivity *)
ϕPrimeRentingFun[η_,q_,ϕ_,d_]:=With[{ψ=Ψ[η,q]},-d+((ϕ (1-(q^2)-2 q η+2 θ (ab-δ+μ)+2 η (1+θ (r+μ))+2 θ ψ (a-ab+η ((σ2)^2) ψ)))/(θ (((σ1)^2)+(η^2) ((σ2)^2)) (ψ^2)))];

ϕPrimeNotRentingFun[η_,q_,ϕ_,d_]:=-d+(((1-q (q+2 η)+2 θ (a-δ+μ)+2 η (1+θ (r+μ+((σ2)^2)))) ϕ)/(θ (((σ1)^2)+(η^2) ((σ2)^2))));
(*Now define two Fokker-Planck equations for subsequent numerical solution, the first when there is no liquidation so 0= -d *)
ϕ0PrimeRentingEq={ϕ0'[η]==ϕPrimeRentingFun[η,q[η],ϕ0[η],0]};

(* the second when there is liquidation, and we choose an arbitrary rate of drift -1 = -d *)
ϕ1PrimeRentingEq={ϕ1'[η]==ϕPrimeRentingFun[η,q[η],ϕ1[η],-1]};

(*No Rent Case when 0= -d *)
ϕ0PrimeNotRentingEq={ϕ0'[η]==ϕPrimeNotRentingFun[η,q[η],ϕ0[η],0]};

(*No Rent Case when -1 = -d *)
ϕ1PrimeNotRentingEq={ϕ1'[η]==ϕPrimeNotRentingFun[η,q[η],ϕ1[η],-1]};

(* The upper boundary of our solver in (η, q) space is the dividend paying curve (DPC). The following help us "orient" ourselves with respect to it: We compute the DPC itself as a solution to an implicit equation, and we solve the η  *)

(* The DPC lies on the curve where q derivative vanishes, yielding the following implicit equation *)
DPCImplicitEq=Simplify[(((σ1)^2)+(η^2) ((σ2)^2))qPrimeNotRentingFun[η,q]/(q+η)==0];

(* Solve the value of ηStar given a value of q. *)
ηStarOfq=η/.First@Solve[DPCImplicitEq,η];

(* Extremal point of the upper boundary curve *)
DPCTip={ηStarOfq,q}/.First@Solve[(∂_q)ηStarOfq==0,q];

(* Solve the value of q on the DPC corresponding to given η *)
qStarsOfη=q/.Solve[DPCImplicitEq,q]; 

(* The minimum η where the DPC still exists (the DPC is a parabola opening to the right on (η,q) plane, so there is a minimum η for which it still exists. *)
DPCηMinimum=η/.First@Solve[Discriminant[qPrimeNotRentingFun[η,q]/(q+η),q]==0,η]; 

(* Value of q at η minimum, hitting the DPC above this point implies we are on the upper branch (which we generally do not want here) *)
qDPCηMinimum=q/.First@Solve[DPCImplicitEq/.η->DPCηMinimum,q]; 

PrintDefValue["DPCηMinimum"]
PrintDefValue["qDPCηMinimum"]

ComputeηSingular[usePars_]:=Module[{ηSing},
ηSing=ηSingular/.usePars;
If[Im@ηSing≠0,Throw[{"ComputeηSingular: Initial η does not exist -- check parameters"},"QSolver"]];
ηSing
]

(*The following two functions compute "initial" values of psi, q' , for the case fo rental. These are based on the assymptotic approximations using ηsingular as the intial arbitrary η 
This is based on equations (35) on page47 of the paper and the unnumbered equation for psiprime on page 48.
As it turns out these "initial" values are also the final ones, provided only that the firm does not choose recapitalisation
  *)
(* MOO *)
ComputeqPrimeInit[usePars_]:=Module[{ηInit,c=(((a-ab))^2)/(2 (((σ1)^2)+(η^2) ((σ2)^2)))},
ηInit=ComputeηSingular[usePars];
((-(ρ-r-c)+(Sqrt[(((ρ-r-c))^2)+((4c)/θ)(1+θ(ρ+μ)+η)]))/(2(θ^(-1))(1+θ(ρ+μ)+η)))/.η->ηInit/.usePars
];  
qPrimeInitDef=ComputeqPrimeInit[defPars];
PrintDefValue["qPrimeInitDef"]

ComputepsiPrimeInit[usePars_]:=Module[{qp,η},
η=ComputeηSingular[usePars];
qp=ComputeqPrimeInit[usePars];
(2/(a-ab))(ρ-r+(1/θ)(1+θ(ρ+μ)+η)qp)/.usePars
] 
ψPrimeInitDef=ComputepsiPrimeInit[defPars];
PrintDefValue["ψPrimeInitDef"];

(* We then have functions for calculating the asymtptotic (power law) approximations, for f() and for W, in the case of rental and also corresponding remainders *) 
ComputefPowerLawExp[usePars_]:=Module[{ψp,ηb,qp},
ηb=ComputeηSingular[usePars];
ψp=ComputepsiPrimeInit[usePars];
qp=ComputeqPrimeInit[usePars];
(((a-ab)ψp+r+μ-(1/(2θ))qp(ηb+1)+(1/θ)(ηb+1)(1+qp/2))/((1/2)(((σ1)^2)+((ηb)^2)((σ2)^2))((ψp)^2)))-2/.usePars
]

ComputeWPowerLawExp[paramRule_]:=1/(1+ComputeqPrimeInit[paramRule])
WPowerLawExpDef=ComputeWPowerLawExp[defPars];
PrintDefValue["WPowerLawExpDef"]

(* Arbitrary max. value for η *)
ηMax=5.0; 

(* This is the most important auxilliary function: It solves the problem for a given initial value of q, generating all the functions q(η), f(η), F(η), and W(η). *)
(* Error behaviour is as follows: If integrator (NDSolve) fails between (OverBar[η]) and (OverTilde[η]), then throw error. If on the other hand error occurs between (OverTilde[η]) and (SuperStar[η]), or the upper boundary condition is not satisfied, or if the upper boundary occurs on the upper branch, then accept the solution but return a status value indicating what happened. *)
statusNoBoundary=1;
statusNoBoundaryNegDrift=2;
statusFailBeforeBoundary=3;
statusUpperBranch=4;
statusOK=0;
AuxSolveqWfF[qI_?NumericQ,paramRule_]:=Module[{
ηI,ηCα,ηCβ,ηC,ηT,ηS,ηR,α,β,qCn,ϵ,Ψn,Σn,ΣNRn,
haveηT=False,haveηS=False,status=statusOK,
qRSEqn,qREqn,qNREqn,ϕ0REqn,ϕ1REqn,ϕ0NREqn,ϕ1NREqn,
qRn,ϕ0Rn,ϕ1Rn,F0Rn,F1Rn,WRn,
qNRn,ϕ0NRn,ϕ1NRn,F0NRn,F1NRn,WNRn,Wn,
a0,a1,
WNorm,FS,
qIn,ϕIn,FIn,WIn,
qn,Fn,fn,
qCrit,ran
},
Σn[η_,q_]=halfσηSquare[η,q]/.paramRule;
ΣNRn[η_,q_]=halfσηSquareNotRenting[η,q]/.paramRule;
Ψn[η_,q_]=Ψ[η,q]/.paramRule;
ϵ=ϵCut/.paramRule;
qRSEqn=qPrimeRentingNearSingularEq/.paramRule;
qREqn=qPrimeRentingEq/.paramRule;
ϕ0REqn=ϕ0PrimeRentingEq/.paramRule;
ϕ1REqn=ϕ1PrimeRentingEq/.paramRule;
qNREqn=qPrimeNotRentingEq/.paramRule;
ϕ0NREqn=ϕ0PrimeNotRentingEq/.paramRule;
ϕ1NREqn=ϕ1PrimeNotRentingEq/.paramRule;
qCrit=qDPCηMinimum/.paramRule;
(* Initial values *)
ηR=ηI=ComputeηSingular[paramRule];
ηC=ηI;
qIn=qI;
ϕIn=1.0;
FIn=0.0;
WIn=1.0;
(* Find cross-over η, point where to switch between approximate analytic and numerical solutions *)
If[qI==-ηI,
(******* IF INITIAL POINT IS THE SINGULAR ONE *********)

(* Need the analytic formulas *)

α=ComputeWPowerLawExp[paramRule];
β=ComputefPowerLawExp[paramRule];
ηCα=ϵ Abs[α]+ηI;
ηCβ=ϵ Abs[β]+ηI;
ηC=Min[ηCα,ηCβ];
(* Extra checks *)
If[ηC-ηI==0.0,Throw[{"AuxSolveqWfF: To within numerical precision, initial η and analytic cross over η are the same. Because of this, cannot construct matching between analytic and numerical solutions."},"QSolver"]];
(* Solve q upto cross-over and verify that we are still in the renting range *)
qCn=q/.First@Quiet@Check[NDSolve[qRSEqn~Join~{q[ηI]==qI},q,{η,ηI,ηC}],Throw[{"AuxSolveqWfF: ODE solver failed while attempting solve q near a singular initial point, ie. point such that ψ=0, q=-η. This error is not expected for typical parameter values. Check parameters, or contact authors."},"QSolver"]];
If[Ψn[ηC,qCn[ηC]]>1.0,Throw[{"AuxSolveqWfF: Solver attempted to use an analytic formula near the initial point. However, the η range that was used for the formula was too large, making the end point lie in a region where ψ < 1 constraint binds. Try lowering the parameter ϵCut, or slightly altering the parameter values. This error is expected to occur only for very limited sets of parameters. If error persists, please contact authors."},"QSolver"]];
qIn=qCn[ηC];
ϕIn=Σn[ηI,qIn](((ηC-ηI))^β);
FIn=(1/(β+1))(((ηC-ηI))^(β+1));
WIn=(((ηC-ηI))^α);
,(* else *)
(********* IF INITIAL POINT NOT AT THE SINGULAR INTERSECTION *********)
QCN = qCn=Function[η,qI]; (* q evaluated at the critical value for η*)

];
(* Now solve until ηTilde *)
If[Ψn[ηC,qIn]<1.0,
(******* INITIAL POINT IN R REGIME ********)
Quiet@Check[
{qRn,ϕ0Rn,ϕ1Rn,F0Rn,F1Rn,WRn}={q,ϕ0,ϕ1,F0,F1,W}/.First@NDSolve[
Flatten[{
qREqn,{q[ηC]==qIn},
ϕ0REqn,{ϕ0[ηC]==ϕIn},
ϕ1REqn,{ϕ1[ηC]==ϕIn},
{F0'[η]==1/Σn[η,q[η]]ϕ0[η],F0[ηC]==FIn},
{F1'[η]==1/Σn[η,q[η]]ϕ1[η],F1[ηC]==FIn},
{W'[η]==W[η]/(η+q[η]),W[ηC]==WIn}
}],{q,ϕ0,ϕ1,F0,F1,W},{η,ηC,ηMax},
Method->{
"EventLocator",
"Event"->(Ψn[η,q[η]]-1),
"EventAction":>Throw[haveηT=True;ηT=η, "StopIntegration"]
}
],
(* If NDSolve posts messages (errors, or whatever) abort, and throw an error message *)
(* ran is the η range where the solver was sucessful. *)
ran=qRn[[1,1]];
Throw[{"AuxSolveqWfF: ODE solver failed when solving for q, W, and the ergodic density in an η range where 0<ψ<1. This is expected, if parameters are such that the upper boundary curve is never hit. Please check parameters, and contact authors if needed."(*<>"[Debugging information: Fail η="<>ToString[ran[[2]]]<>", (q,ϕ0,ϕ1,F0,F1,W)="<>ToString[InputForm/@{qRn[ran[[2]]],ϕ0Rn[ran[[2]]],ϕ1Rn[ran[[2]]],F0Rn[ran[[2]]],F1Rn[ran[[2]]],WRn[ran[[2]]]}]<>"]"*)},"QSolver"]];

(*If[Ψn[ηC,qIn]<1.0,Print["Ψn[η,q] < 1 : ",Ψn[ηC,qIn]," and χ is ",χ],Print["Ψn[η,q] NOT < 1 : ",Ψn[ηC,qIn]," and χ is ",χ]];*)

If[!haveηT,Throw[{"AuxSolveqWfF: Solver failed to locate the point where ψ becomes 1. This should always happen before reaching η*, the range end point. A possible reason is that that parameters make the η range unexpectedly long. Try modifying the ηMax parameter (this is a built-in variable, and requires changing this notebook itself, contact authors for help if needed)."},"QSolver"]];
,(* else *)

(********** INITIAL POINT IN NR ********)
ηT=ηI;
qRn=Function[η,qIn];WRn=Function[η,WIn];
ϕ0Rn=Function[η,ϕIn];ϕ1Rn=Function[η,ϕIn];
F0Rn=Function[η,FIn];F1Rn=Function[η,FIn];
];
(* Then continue upto ηStar *)
Quiet@Check[
{qNRn,ϕ0NRn,ϕ1NRn,F0NRn,F1NRn,WNRn}={q,ϕ0,ϕ1,F0,F1,W}/.First@NDSolve[
Flatten[{
qNREqn,{q[ηT]==qRn[ηT]},
ϕ0NREqn,{ϕ0[ηT]==ϕ0Rn[ηT]},
ϕ1NREqn,{ϕ1[ηT]==ϕ1Rn[ηT]},
{F0'[η]==1/ΣNRn[η,q[η]]ϕ0[η],F0[ηT]==F0Rn[ηT]},
{F1'[η]==1/ΣNRn[η,q[η]]ϕ1[η],F1[ηT]==F1Rn[ηT]},
{W'[η]==W[η]/(η+q[η]),W[ηT]==WRn[ηT]}
}],{q,ϕ0,ϕ1,F0,F1,W},{η,ηT,ηMax},
Method->{
"EventLocator",
"Event"->q'[η],
"EventAction":>Throw[haveηS=True;ηS=η, "StopIntegration"]
}
],
(* This is evaluated if NDSolve generated messages. In this case, don't throw an error, but do mark the status flag as !statusOK. *)
status=statusFailBeforeBoundary;
haveηS=True;
ηS=qNRn[[1,1,2]];
(* OLD handler, this simply threw an error *)
(* ran=qNRn[[1,1]];
Throw[{"AuxSolveqWfF: ODE solver failed when solving for q and the ergodic density in an η range where ψ=1. This error is unexpected as no singularities are present. This error may occur if given parameters do not have numerical values, or are otherwise too extreme. Please check parameters, and contact authors if needed."<>"[Debugging information: Fail η="<>ToString[ran[[2]]]<>", (q,ϕ0,ϕ1,F0,F1,W)="<>ToString[InputForm/@{qNRn[ran[[2]]],ϕ0NRn[ran[[2]]],ϕ1NRn[ran[[2]]],F0NRn[ran[[2]]],F1NRn[ran[[2]]],WNRn[ran[[2]]]}]<>"Init vals: η="<>ToString[ran[[1]]]<>", state="<>ToString[InputForm/@{qNRn[ran[[1]]],ϕ0NRn[ran[[1]]],ϕ1NRn[ran[[1]]],F0NRn[ran[[1]]],F1NRn[ran[[1]]],WNRn[ran[[1]]]}]<>"]"},"QSolver"]*)
];
If[!haveηS,
(* Boundary not reached, no other errors. Set status to reflect this and the region where the integration ended (negative vs. positive drift) *)
haveηS=True;ηS=ηMax;
If[qNRn[ηMax]>qCrit
,status=statusNoBoundaryNegDrift
,status=statusNoBoundary
];
,
If[qNRn[ηS]>qCrit&&status==statusOK,status=statusUpperBranch];
];
(* Normalize W *)
WNorm=((WNRn'[ηS])^(-1));
WNORM = WNorm;(*EXPOSE THE VARIABLE *)
Wn=Function[η,Evaluate[WNorm({{\[Piecewise], ({{(((η-ηI))^α), ηI≤η<ηC}, {WRn[η], ηC≤η<ηT}, {WNRn[η], ηT≤η≤ηS}})}})]];
WN = Wn; (*EXPOSE THE VARIABLE *)
(* Construct the normalised pdf and cdf *)
If[qI==-ηI,
(*** Reflecting state ***)
a0=1/F0NRn[ηS];a1=0.0;
REFLECTING = True;
ABSORBING =False;
fn=If[β<-1.0,DiracDelta[#-ηI]&,Function[η,Evaluate[({{\[Piecewise], ({{(a0(((η-ηI))^β)), ηI≤η<ηC}, {1/Σn[η,qRn[η]](a0 ϕ0Rn[η]), ηC≤η<ηT}, {1/ΣNRn[η,qNRn[η]](a0 ϕ0NRn[η]), ηT≤η≤ηS}})}})]]];
Fn=If[β<-1.0,UnitStep[#-ηI]&,Function[η,Evaluate[a0({{\[Piecewise], ({{1/(β+1)(((η-ηI))^(β+1)), ηI≤η<ηC}, {F0Rn[η], ηC≤η<ηT}, {F0NRn[η], ηT≤η≤ηS}})}})]]]
,(* else *)
(*** Absorbing state ***)
REFLECTING = False;
ABSORBING =True;
a0=ϕ1Rn[ηI]/(-F1NRn[ηS] ϕ0Rn[ηI]+F0NRn[ηS] ϕ1Rn[ηI]);
a1=ϕ0Rn[ηI]/(F1NRn[ηS] ϕ0Rn[ηI]-F0NRn[ηS] ϕ1Rn[ηI]);
fn=Function[η,Evaluate[({{\[Piecewise], ({{1/Σn[η,qRn[η]](a0 ϕ0Rn[η]+a1 ϕ1Rn[η]), ηI≤η<ηT}, {1/ΣNRn[η,qNRn[η]](a0 ϕ0NRn[η]+a1 ϕ1NRn[η]), ηT≤η≤ηS}})}})]];
Fn=Function[η,Evaluate[({{\[Piecewise], ({{a0 F0Rn[η]+a1 F1Rn[η], ηI≤η<ηT}, {a0 F0NRn[η]+a1 F1NRn[η], ηT≤η≤ηS}})}})]]
];
qn=Function[η,Evaluate@({{\[Piecewise], ({{qCn[η], ηI≤η<ηC}, {qRn[η], ηC≤η<ηT}, {qNRn[η], ηT≤η≤ηS}})}})];

ETAR =ηR;
ETACRITICAL =ηC;
ETATILDER = ηT;
ETASTAR =ηS;
QN = qn;
fN =fn;
FN = Fn;
WN = Wn;


(* Return the critical η values and piecewise solutions *)
{
status,
ηR,ηC,ηT,ηS,
qn,fn,Fn,Wn
}
]

(* Back-substitution of solution of Fokker-Planck. Takes a solution of f, fp = f' of the FP equation, and evaluates the differential equation. If this yields zero (or constant), we are good (the solution satisfies the equation) *)
(* NOT IMPLENTED YET FOR μ ≠ 0 *)
(* 
CompfPrimeRemainder[η_,q_,qp_,f_,fp_,usePars_]:=Module[{ψn,dψn,σn,dσn,μn},
ψn=Ψ[η,q]/.usePars;
dψn=If[0≤ψn<1,-(2 (r+1/θ(-1+q) qp-ρ (1+qp)))/(a-ab),0]/.usePars;
ψn=Clip[ψn,{0,1}];
σn=√(((σ1)^2)+η ((σ2)^2))ψn/.usePars;
dσn=(((σ2)^2) ψn)/(2 √(((σ1)^2)+η ((σ2)^2)))+√(((σ1)^2)+η ((σ2)^2))dψn/.usePars;
μn=ab+(a-ab)ψn-δ+r η-(θ^(-1))(1+η)(q-1)-(1/2)(θ^(-1))(((q-1))^2)+η ((σ2)^2)((ψn)^2)/.usePars;
μn f-dσn σn f-(1/2)((σn)^2)fp
];
*)

(*Now the "recapitalizationGain" based on eqn (12) page 18 of the paper, if this is positive the firm recapitalises rather than liquidates *)
(* This uses AuxSolveqWfF. If error is caught, it is simply thrown back, with an additional error messages. If boundary conditions are not satisfied, accept the solution but also return the status code. *)
Clear[ComputeRecapitalizationGain];
ComputeRecapitalizationGain[q0_?NumericQ,usePars_]:=Module[{status,ηIx,ηRx,ηCx,ηTx,ηSx,qnx,fnx,Fnx,Wnx},
ηIx=ComputeηSingular[usePars];
Catch[
{status,ηRx,ηCx,ηTx,ηSx,qnx,fnx,Fnx,Wnx}=AuxSolveqWfF[q0,usePars];
,"QSolver"
,Throw[{"ComputeRecapitalizationGain: Failed computing the gain from recapitalization, as the solver method is reporting an error."}~Join~#1,"QSolver"]&
];
RECAPGAIN = Wnx[ηSx]-Wnx[ηIx]-(ηSx-ηIx)-(χ/.usePars);
CHI = χ/.usePars;
{status,Wnx[ηSx]-Wnx[ηIx]-(ηSx-ηIx)-(χ/.usePars)}
]

(* This function computes an upper bound for q at initial η. Used by function ComputeLowerBoundaryq to setup a search interval for initial q. The upper bound is either (i) the dividend paying curve, if the dpc intersects the initial η value, or (ii) the q value for which the branches of the dpc meet (as there is no point in searching higher values of q: they are bound to hit the dpc on the upper branch). *)
ComputeqUpperBound[usePars_]:=Module[{qmax,η0,ηdpc,qdpc},
ETAZERO =η0=ComputeηSingular[usePars];
ETADPC =ηdpc=DPCηMinimum/.usePars;
If[η0<ηdpc,
{qDPCηMinimum/.usePars,False},
{Min[qStarsOfη/.usePars/.η->η0],True}
]
] 
qUpperBoundDef=ComputeqUpperBound[defPars];
PrintDefValue["qUpperBoundDef"];


(* Compute W and q on the upper boundary for given initial q, q0. Auxilliary function used for computing the lower bound q. *)
ComputeLocalWqStar[q0_?NumericQ,usePars_]:=Module[{status,ηIx,ηCx,ηTx,ηSx,qnx,fnx,Fnx,Wnx},
{status,ηIx,ηCx,ηTx,ηSx,qnx,fnx,Fnx,Wnx}=Catch[AuxSolveqWfF[q0,usePars],"QSolver",Throw[{"ComputeLocalWqStar: Attempted to compute value function and q for some q initial condition, however, solver failed"}~Join~#1,"QSolver"]&];
{WNXSTAR,QNXSTAR} = {Wnx[ηIx],qnx[ηSx]}  ;
{status,Wnx[ηIx],qnx[ηSx]}
];

(* the following function is used within "solveverything" to search for the initial value of q  that  satisfies the matching condition (equatino (12) on page 18 of the Nov 2012 paper if this is satsified, or if not liquidation .

it is relatively complicated because of the possibility of "lower branch" v. "upper branch"

the basic calculation however is very simple, numericall search on q(etabar) to find the value that maximises W(etabar), this is when G gain fro increasing q(etabar) by a small amount is -ve
   but Gprev gain from reducing q(etabar) by a small amount is +ve)
 *)
(* This function scans a range of q values at the initial η, ranging from -η to an upper boundary specified in function ComputeqUpperBound. A special case: If the solver fails at the lowest possible q value, for whatever reason *)
ComputeLowerBoundaryq[usePars_]:=Module[{status,η0,ηSx,q0x,q1x,nqq=qBarSearchNDiv/.usePars,j,n,dqq,qfound,qstar,qmax,G,Gprev,qbest,Wbar,Wbest=0.0,dpcbinds,χn=χ/.usePars,recap=False,q0xShift,haveroot,fun},
qmax=qDPCηMinimum/.usePars;
η0=ComputeηSingular[usePars];
(* Compute the bounds for the initial/maxdebt q variable search *)
QZEROX = q0x=-η0;
{q1x,dpcbinds}=ComputeqUpperBound[usePars];
(* Search for roots *)
dqq=(q1x-q0x)/nqq;
(* There's a slight danger of failure here: If a root is bracketed on the first iteration, the root finder might take us too close to the singular point, leading to the integrator failing. What we do is we quietly accept failure of the root finder if the failure occurs on the first step of bracketing. In this case we have succesfully evaluated the recapitalisation gain at both ends of a small bracketing interval, the beginning being the singular point, and found that a root is somewhere in this interval. If the root finder/polisher fails in this case, we take it to mean that the root is so close to the singular point that we cannot evaluate the gain. It should be then reasonable to take q = -(Subscript[η, singular]) as the initial point, and suppose liquidation. *)

{status,Gprev}=Catch[ComputeRecapitalizationGain[q0x,usePars],"QSolver",Throw[{"ComputeLowerBoundaryq: Could not compute the lower boundary for q, as finding the gain from recapitalization failed at lowest q possible"}~Join~#1,"QSolver"]&];

(* If the solver fails at the lowest possible value of q by failing to hit the upper boundary, then the boundary is not expected to be hit for any lower value either. We take the -η as the initial q, and return a status indicating the abnormal boundary. *)
If[status≠statusOK,
Return[{status,q0x}]
];

(* Wrap the ComputeRec...Gain into a function that just throws errors on bad status cases. This is more digestible to FindRoot. *)
fun[qq_?NumericQ]:=Module[{stat,GG},{stat,GG}=ComputeRecapitalizationGain[qq,usePars];If[stat≠statusOK,Throw["Fail","QSolver"]];GG];

Catch[
For[j=1,j≤nqq,j++,
If[j==nqq&&dpcbinds
,status=statusOK;G=-χn;
,{status,G}=ComputeRecapitalizationGain[q0x+j dqq,usePars];
];
(* If this q gives abnormal boundary behaviour, then quit search loop. *)
If[status≠statusOK,Throw["Fail","QSolver"]];
(*Print["G=",G,"Gprev=",Gprev];*)
If[G Gprev≤0,
haveroot=False;
(*Print["Finding root, q in (",q0x+(j-1)dqq,",",q0x+j dqq,")"];*)
Catch[
Check[
qfound=qq/.FindRoot[fun[qq]==0,{qq,q0x+(j-1)dqq,q0x+j dqq}];
,Throw[{"Fail"},"QSolver"] (* This throw is caught before it gets to the user, no need for a long error message yet. *)
];
haveroot=True;
,"QSolver"
];
If[!haveroot&&j>1,Throw[{"ComputeLowerBoundaryq: Failed to find initial values satisfying the racapitalization matching condition. This error is unexpected as initial, rough bracketing indicates that the solution exists, but root polishing still failed. Contact authors for assistance."},"QSolver"]];
(* We have found a q value satisfying matching. 
* Check that this solution does not yield an investment policy that hits DPC on its upper branch. *)
If[haveroot,
{status,Wbar,qstar}=ComputeLocalWqStar[qfound,usePars];
If[status==statusOK&&Wbar>Wbest&&qstar<qmax,
recap=True;
qbest=qfound;
Wbest=Wbar;
];
];
];
Gprev=G;
];
,"QSolver"
];
If[recap,
(* Recapitalize. Status is verified to be OK, so return: *)
Return[{statusOK,qbest}],
(* Before accepting liquidation, just find the upper boundary behaviour, and pass it to caller: *)
{status,Wbar,qstar}=Catch[ComputeLocalWqStar[q0x,usePars],"QSolver",Throw[{"ComputeLowerBoundaryq: Recapitalization found not possible, but could not compute value function and q to check validity of the solution." }~Join~#1]&];
Return[{status,q0x}];
];
]
lowerBoundaryqDef=ComputeLowerBoundaryq[defPars];
PrintDefValue["lowerBoundaryqDef"];

SolveEverything[paramRule_]:=Module[{status,qIx,ηIx,ηCx,ηTx,ηSx,qnx,fnx,Fnx,Wnx},
Catch[
{status,qIx}=ComputeLowerBoundaryq[paramRule];
,"QSolver",Throw[{"SolveEverything: Could not determine the lower boundary q value."}~Join~#1,"QSolver"]&
];
{status,ηIx,ηCx,ηTx,ηSx,qnx,fnx,Fnx,Wnx}=AuxSolveqWfF[qIx,paramRule]
];

Model without rental

(************ Set up of solution routines for Model without rent ************)
(* Equations (some repetition with the paramterisation cell above ) ) *)
qPrimeNoRentFun[η_,q_]=(a-δ-(ρ-r)η-ρ q+(1/2)(θ^(-1))(((q-1))^2))/((1/2)((σ1)^2))(q+η);
qPrimeNoRentEq={q'[η]==1/((1/2)((σ1)^2))(a-δ-(ρ-r)η-ρ q[η]+(1/2)(θ^(-1))(((q[η]-1))^2))(q[η]+η)};

WPrimeNoRentEq={W'[η]==If[q[η]+η==0.0,1.0,W[η]/(q[η]+η)]};

ηDriftNoRentFun[η_,q_]:=a-δ+r η-(1+η)(θ^(-1))(q-1)-(1/2)(θ^(-1))(((q-1))^2);
ηDiffusionNoRentFun[η_,q_]:=(1/2)((σ1)^2);
fDriftNoRentFun[η_,q_]:=ηDriftNoRentFun[η,q]/ηDiffusionNoRentFun[η,q];

fPrimeNoRentEq={f'[η]==fDriftNoRentFun[η,q[η]]f[η]};

fPrimeInhomogNoRentEq={fp'[η]==fDriftNoRentFun[η,q[η]]fp[η]+1}; 


(* Numerical solver *)
AuxSolveqWfFNoRent[qb_,testParameters_]:=Module[{qInit,ηInit,ηStar,ηR,qSol,fSol,FSol,fpSol,FpSol,fRSol,FRSol,WSol,equation,fReflSol,FReflSol,fAbsSol,FAbsSol,WSolN,normR,a0,a1,Wnorm,qCrit,haveηStar=False,χn,recapitalize,status=statusOK},
χn=χ/.testParameters;
(* Initial values *)
ηR=ηInit=ComputeηSingular[testParameters];
qInit=qb;
qCrit=qDPCηMinimum/.testParameters;
RECAPNR =recapitalize=!(qInit==-ηInit);
(* Construct the initial value problem *)
equation=Flatten@{
qPrimeNoRentEq,q[ηInit]==qInit,
fPrimeNoRentEq,f[ηInit]==1,{F'[η]==f[η],F[ηInit]==0},
WPrimeNoRentEq,W[ηInit]==If[recapitalize,1,0],
fPrimeInhomogNoRentEq,fp[ηInit]==1,{Fp'[η]==fp[η],Fp[ηInit]==0}
}/.testParameters;
(* Solve it *)
Quiet@Check[
{qSol,fSol,FSol,fpSol,FpSol,WSol}={q,f,F,fp,Fp,W}/.First@NDSolve[equation,{q,f,F,fp,Fp,W},{η,ηInit,ηMax},
Method->{
"EventLocator",
"Event"->q'[η],
"EventAction":>Throw[haveηStar=True;ηStar=η, "StopIntegration"]
}
,MaxSteps->(10^6)
],
(* This is evaluated if NDSolve generated messages. In this case, don't throw an error, but do mark the status flag as !statusOK. *)
status=statusFailBeforeBoundary;
haveηStar=True;
ηStar=qSol[[1,1,2]];
(* *)
] (* Check *);
(* Success, but is ηStar found? *)
If[!haveηStar,
(* Boundary not reached, no other errors. Set status to reflect this and the region where the integration ended (negative vs. positive drift) *)
haveηStar=True;ηStar=ηMax;
If[qSol[ηMax]>qCrit
,status=statusNoBoundaryNegDrift
,status=statusNoBoundary
];
,
If[qSol[ηStar]>qCrit&&status==statusOK,status=statusUpperBranch];
];
(* Normalize W so that W'(ηStar)==1 *)
WNORMNR = Wnorm=WSol'[ηStar];
WSolN=Function[η,WSol[η]/Wnorm];
(* Construct normalized densities *)
(* Absorbing boundary  *)
a0=fpSol[ηR]/(FSol[ηStar] fpSol[ηR]-fSol[ηR] FpSol[ηStar]);
a1=fSol[ηR]/(-FSol[ηStar] fpSol[ηR]+fSol[ηR] FpSol[ηStar]);
fRSol=Function[η,a0 fSol[η]+a1 fpSol[η]];
FRSol=Function[η,a0 FSol[η]+a1 FpSol[η]];
{ETARNR,ETASTARNR,QSOL,WSOLN,fRSOL,FRSOL}={ηR,ηStar,qSol,WSolN,fRSol,FRSol};
{status,ηR,ηStar,qSol,fRSol,FRSol,WSolN}
];

ComputeRecapitalizationGainNoRent[q0_?NumericQ,usePars_]:=Module[{status,ηIx,ηRx,ηCx,ηTx,ηSx,qnx,fnx,Fnx,Wnx},
ηIx=ComputeηSingular[usePars];
Catch[
{status,ηRx,ηSx,qnx,fnx,Fnx,Wnx}=AuxSolveqWfFNoRent[q0,usePars];
,"QSolver"
,Throw[{"ComputeRecapitalizationGainNoRent: Failed computing the gain from recapitalization, as the solver method is reporting an error."}~Join~#1,"QSolver"]&
];
RECAPGAINNR = Wnx[ηSx]-Wnx[ηIx]-(ηSx-ηIx)-(χ/.usePars);
CHI = χ/.usePars;
{status,Wnx[ηSx]-Wnx[ηIx]-(ηSx-ηIx)-(χ/.usePars)}
]

ComputeLocalWqStarNoRent[q0_?NumericQ,usePars_]:=Module[{status,ηIx,ηCx,ηTx,ηSx,qnx,fnx,Fnx,Wnx},
{status,ηIx,ηSx,qnx,fnx,Fnx,Wnx}=Catch[
AuxSolveqWfFNoRent[q0,usePars]
,"QSolver"
,Throw[{"ComputeLocalWqStarNoRent: Attempted to compute value function and q for some q initial condition, however, solver failed"}~Join~#1,"QSolver"]&
];
{WNXSTARNR,QNXSTARNR}={Wnx[ηIx],qnx[ηSx]};
{status,Wnx[ηIx],qnx[ηSx]}
];

ComputeLowerBoundaryqNoRent[usePars_]:=Module[{status,η0,ηSx,q0x,q1x,nqq=qBarSearchNDiv/.usePars,j,n,dqq,qfound,qstar,qmax,G,Gprev,qbest,Wbar,Wbest=0.0,dpcbinds,χn=χ/.usePars,recap=False,q0xShift,haveroot,fun},
qmax=qDPCηMinimum/.usePars;
η0=ComputeηSingular[usePars];
(* Compute the bounds for the initial/maxdebt q variable search *)
QZEROX = q0x=-η0;
{q1x,dpcbinds}=ComputeqUpperBound[usePars];
(* Search for roots *)
dqq=(q1x-q0x)/nqq;
(* There's a slight danger of failure here: If a root is bracketed on the first iteration, the root finder might take us too close to the singular point, leading to the integrator failing. What we do is we quietly accept failure of the root finder if the failure occurs on the first step of bracketing. In this case we have succesfully evaluated the recapitalisation gain at both ends of a small bracketing interval, the beginning being the singular point, and found that a root is somewhere in this interval. If the root finder/polisher fails in this case, we take it to mean that the root is so close to the singular point that we cannot evaluate the gain. It should be then reasonable to take q = -(Subscript[η, singular]) as the initial point, and suppose liquidation. *)

{status,Gprev}=Catch[ComputeRecapitalizationGainNoRent[q0x,usePars],"QSolver",Throw[{"ComputeLowerBoundaryqNoRent: Could not compute the lower boundary for q, as finding the gain from recapitalization failed at lowest q possible"}~Join~#1,"QSolver"]&];

(* If the solver fails at the lowest possible value of q by failing to hit the upper boundary, then the boundary is not expected to be hit for any lower value either. We take the -η as the initial q, and return a status indicating the abnormal boundary. *)
If[status≠statusOK,Return[{status,q0x}]];

(* Wrap the ComputeRec...Gain into a function that just throws errors on bad status cases. This is more digestible to FindRoot. *)
fun[qq_?NumericQ]:=Module[{stat,GG},{stat,GG}=ComputeRecapitalizationGainNoRent[qq,usePars];If[stat≠statusOK,Throw["Fail","QSolver"]];GG];

Catch[
For[j=1,j≤nqq,j++,
If[j==nqq&&dpcbinds
,status=statusOK;G=-χn;
,{status,G}=ComputeRecapitalizationGainNoRent[q0x+j dqq,usePars];
];
(* If this q gives abnormal boundary behaviour, then quit search loop. *)
If[status≠statusOK,Throw["Fail","QSolver"]];
(*Print["G=",G,"Gprev=",Gprev];*)
If[G Gprev≤0,
haveroot=False;
(*Print["Finding root, q in (",q0x+(j-1)dqq,",",q0x+j dqq,")"];*)
Catch[
Check[
qfound=qq/.FindRoot[fun[qq]==0,{qq,q0x+(j-1)dqq,q0x+j dqq}];
,Throw[{"Fail"},"QSolver"] (* This throw is caught before it gets to the user, no need for a long error message yet. *)
];
haveroot=True;
,"QSolver"
];
If[!haveroot&&j>1,Throw[{"ComputeLowerBoundaryqNoRent: Failed to find initial values satisfying the racapitalization matching condition. This error is unexpected as initial, rough bracketing indicates that the solution exists, but root polishing still failed. Contact authors for assistance."},"QSolver"]];
(* We have found a q value satisfying matching. 
* Check that this solution does not yield an investment policy that hits DPC on its upper branch. *)
If[haveroot,
{status,Wbar,qstar}=ComputeLocalWqStarNoRent[qfound,usePars];
If[status==statusOK&&Wbar>Wbest&&qstar<qmax,
recap=True;
qbest=qfound;
Wbest=Wbar;
];
];
];
Gprev=G;
];
,"QSolver"
];
If[recap,
(* Recapitalize. Status is verified to be OK, so return: *)
Return[{statusOK,qbest}],
(* Before accepting liquidation, just find the upper boundary behaviour, and pass it to caller: *)
{status,Wbar,qstar}=Catch[ComputeLocalWqStarNoRent[q0x,usePars],"QSolver",Throw[{"ComputeLowerBoundaryqNoRent: Recapitalization found not possible, but could not compute value function and q to check validity of the solution." }~Join~#1]&];
Return[{status,q0x}];
];
]  (* We obtain 0.5439919110295031`*)


SolveEverythingNoRent[paramRule_]:=Module[{status,qIx,ηIx,ηSx,qnx,Wnx,fnx,Fnx},
{status,qIx}=ComputeLowerBoundaryqNoRent[paramRule];
{status,ηIx,ηSx,qnx,fnx,Fnx,Wnx}=AuxSolveqWfFNoRent[qIx,paramRule]
]

(*************** END OF SET UP OF ROUTINES FOR MODEL WITHOUT RENT **************)

Routines for analysis of solutions

These functions compute additional quantities of interest from the solutions generated by the solver algorithms of previous section.

(* Further routines for computing countour plots of scaled median  *)
ComputeScaledMedian[usePars_]:=Module[{status,ηIx,ηCx,ηTx,ηSx,qnx,fnx,Fnx,Wnx,ηMedian,β},
Catch[
{status,ηIx,ηCx,ηTx,ηSx,qnx,fnx,Fnx,Wnx}=SolveEverything[usePars];
If[status≠statusOK,Return[I]];
If[ηIx==-qnx[ηIx]
,β=ComputefPowerLawExp[usePars];
ηMedian=If[β<-1.0,ηIx,η/.FindRoot[Fnx[η]==1/2,{η,ηIx,ηSx}]]
,ηMedian=η/.FindRoot[Fnx[η]==1/2,{η,ηIx,ηSx}];
];
SCALEDMEDIAN=(ηMedian-ηIx)/(ηSx-ηIx);
(ηMedian-ηIx)/(ηSx-ηIx)
,
"QSolver",
(I&)
]
]

ComputeScaledMedianNoRent[usePars_]:=Module[{status,ηIx,ηCx,ηTx,ηSx,qnx,fnx,Fnx,Wnx,ηMedian,β},
Catch[
{status,ηIx,ηSx,qnx,fnx,Fnx,Wnx}=SolveEverythingNoRent[usePars];
If[status≠statusOK,Return[I]];
ηMedian=η/.FindRoot[Fnx[η]==1/2,{η,ηIx,ηSx}];
ETAMEDIANNR =ηMedian;
SCALEDMEDIANNR =(ηMedian-ηIx)/(ηSx-ηIx);
(ηMedian-ηIx)/(ηSx-ηIx)
,
"QSolver",
(I&)
]
]

Solve-and-plot -routines

(* Now the SolveAndPlotEverythingRoutines for main graphics*)

NF[x_]:={x,Style[Quiet@NumberForm[x,2],Background->White]};

StatusMsg[status_]:=Switch[status
,statusOK,""
,statusNoBoundary,"Error: Upper boundary not reached. As q is below upper branch of boundary curve, this is unexpected."
,statusNoBoundaryNegDrift,"Error: Upper boundary not reached. Since q above upper branch, likely reason is solutions never hitting the upper boundary curve."
,statusFailBeforeBoundary,"Error: The upper boundary condition could not be satisfied. Likely cause is solutions never hitting the upper boundary curve."
,statusUpperBranch,"Warning: Upper branch solution."
];
SolveAndPlotEverything[usePars_]:=Module[{status,ηIx,ηCx,ηTx,ηSx,qnx,fnx,Fnx,Wnx,Wplot,qplot,ψplot,Fplot,mt,qtick,ηtick,ψn},
{status,ηIx,ηCx,ηTx,ηSx,qnx,fnx,Fnx,Wnx}=SolveEverything[usePars];

ψn[η_]=Ψ[η,qnx[η]]/.usePars;
qplot=Plot[
{qnx[η],qnx[ηIx]+(qnx[ηSx]-qnx[ηIx])Clip[ψn[η],{0,1}]},
{η,ηIx,ηSx}
,Frame->True
,Axes->False
,PlotStyle->{{ColorData[1,1]},{ColorData[1,1],Dashed}}
,FrameTicks->{{{NF@qnx[ηIx],NF@(qnx[ηTx]),NF@qnx[ηSx]},{{qnx[ηIx],"0.0"},{qnx[ηSx],"1.0"}}},{{NF@ηIx,NF@(ηTx),NF@ηSx},None}}
(*,FrameTicks->{{Table[0.2],None},{{NF@ηIx,NF@ηSx},None}}*)
,FrameLabel->{{"q(η)",Style["ψ(η)",Background->White]},{None,None}}
,AspectRatio->1
(*,Epilog->{LightGray,Line[{{ηTx,0},{ηTx,1000}}]}*)
,LabelStyle->12
];
Wplot=Plot[
{Wnx[η]},
{η,ηIx,ηSx}
,Frame->True
,Axes->False
,FrameLabel->{"η","W(η)"} 
,PlotRangeClipping->False
,ImagePadding->{{All,20},{All,All}}
,Epilog->({
Text[Style["(a) ",14],Scaled@{-0.11,1},{0,1}]
(*,Text[Style[StatusMsg[status],14],Scaled@{0.1,0.9},{-1,1}]*)
(*Text[Style["(a)",14],Scaled@{0.1,0.85}]*)
,{LightGray,Line[{{ηTx,0},{ηTx,1000}}]}
,Inset[qplot,Scaled@{0.925,0.125},Scaled@{1,0},Scaled@{0.65,0.65}]
}~Join~{If[status≠statusOK,Inset[Panel@StatusMsg[status],Scaled@{0.05,0.95},{-1,1},{0.4(ηSx-ηIx),Automatic}],{}]})
,LabelStyle->12
];
Fplot=Plot[
{fnx[η],Fnx[η]},
{η,ηIx,ηSx}
,Frame->True
,Axes->False
,FrameLabel->{"η","f(η), F(η)"}
,PlotStyle->{{ColorData[1,1]},{ColorData[1,1],Dashed}}
,LabelStyle->12
,PlotRangeClipping->False
,ImagePadding->{{All,20},{All,All}}
,Epilog->{
Text[Style["(b)",14],Scaled@{-0.11,1},{0,1}]
,{LightGray,Line[{{ηTx,0},{ηTx,1000}}]}
,Text[Style["η = (OverTilde[η])",12],{ηTx,1.0},{0,-1},{0,1}]
(*,Inset[ψplot,Scaled@{0.95,0.15},Scaled@{1,0},Scaled@{0.75,0.75}]*)
}
(*,Epilog->Point[{(ηSx-ηIx)mt+ηIx,1/2}]*)
];
GraphicsGrid[(({{Wplot}, {Fplot}}))]
]



SolveAndPlotEverythingNoRent[usePars_]:=Module[{status,ηIx,ηCx,ηTx,ηSx,qnx,fnx,Fnx,Wnx,Wplot,qplot,Fplot,mt,qtick,ηtick,ψn},
{status,ηIx,ηSx,qnx,fnx,Fnx,Wnx}=SolveEverythingNoRent[usePars];
(*mt=ComputeScaledMedian[usePars];*)
(*qtick=(qnx[ηSx]-qnx[ηSx])/4;*)
qplot=Plot[
{qnx[η]},
{η,ηIx,ηSx}
,Frame->True
,Axes->False
,FrameTicks->{{{NF@qnx[ηIx],NF@((1/2)(qnx[ηIx]+qnx[ηSx])),NF@qnx[ηSx]},None},{{NF@ηIx,NF@((1/2)(ηIx+ηSx)),NF@ηSx},None}}
(*,FrameTicks->{{Table[0.2],None},{{NF@ηIx,NF@ηSx},None}}*)
,FrameLabel->{None,"q(η)"}
,AspectRatio->1
,LabelStyle->12
];
Wplot=Plot[
{Wnx[η]},
{η,ηIx,ηSx}
,Frame->True
,Axes->False
,FrameLabel->{"η","W(η)"}
,PlotRangeClipping->False
,Epilog->{
Text[Style["(a)",14],Scaled@{-0.11,1},{0,1}]
,Inset[qplot,Scaled@{0.925,0.125},Scaled@{1,0},Scaled@{0.65,0.65}]
}
,LabelStyle->12
];
Fplot=Plot[
{fnx[η],Fnx[η]},
{η,ηIx,ηSx}
,Frame->True
,Axes->False
,FrameLabel->{"η","f(η), F(η)"}
,PlotRangeClipping->False
,PlotStyle->{{ColorData[1,1]},{ColorData[1,1],Dashed}}
,LabelStyle->12
,Epilog->{
Text[Style["(b)",14],Scaled@{-0.11,1},{0,1}]
}
(*,Epilog->Point[{(ηSx-ηIx)mt+ηIx,1/2}]*)
];
GraphicsGrid[(({{Wplot}, {Fplot}}))]
]
(*** END OF CELL ****)

Interactive Plots

MakeLabels[xLabel_,{yLabel1_,yLabel2_}]:={xLabel,Row[{yLabel1<>" ",Style["—",ColorData[1,1]],", "<>yLabel2<>" ",Style["—",ColorData[1,2]]}]}
MakeParameterValues[usePars_]:=Module[{},
Text[
Grid[(({{"a", "=", ToString[a/.usePars], "  ", "(OverBar[a])", "=", ToString[ab/.usePars]}, {"δ", "=", ToString[δ/.usePars], "  ", "θ", "=", ToString[θ/.usePars]}, {"r", "=", ToString[r/.usePars], "  ", "ρ", "=", ToString[ρ/.usePars]}, {"(Subscript[σ, 1])", "=", ToString[σ1/.usePars], "  ", "(Subscript[σ, 2])", "=", ToString[σ2/.usePars]}, {"χ", "=", ToString[χ/.usePars], "  ", "μ", "=", ToString[μ/.usePars]}}))
,Alignment->(({{Center, Center, Left, Center, Center, Center, Left}, {Center, Center, Left, Center, Center, Center, Left}, {Center, Center, Left, Center, Center, Center, Left}, {Center, Center, Left, Center, Center, Center, Left}, {Center, Center, Left, Center, Center, Center, Left}}))
,Spacings->{0.2,0.5}
],Scaled@{0.35,0.9},{0,1}]
];
SafetyEpsilon=(10^(-9));
MakePDFCDFPlot[usePars_,f_,F_,{ηInit_,ηStar_}]:=Module[{},
Plot[{f[η],F[η]},{η,ηInit+SafetyEpsilon,ηStar}
,Frame->True
,Axes->False
,FrameLabel->MakeLabels["η",{"f(η)","F(η)"}]
]
]
MakePDFCDFAPlot[usePars_,f_,F_,fA_,FA_,{ηInit_,ηStar_}]:=Module[{},
Show[
Plot[{fA[η],FA[η]},{η,ηInit+SafetyEpsilon,ηStar}
,PlotStyle->{{Gray},{Gray,Dashed}}
,Frame->True
,Axes->False
,FrameLabel->MakeLabels["η",{"f(η)","F(η)"}]
],
Plot[{f[η],F[η]},{η,ηInit+SafetyEpsilon,ηStar}
,Frame->True
,Axes->False
,FrameLabel->MakeLabels["η",{"f(η)","F(η)"}]
]
]
]
MakeDensityPlot[usePars_,f_,{ηInit_,ηStar_}]:=Module[{},
Plot[{f[η]},{η,ηInit+SafetyEpsilon,ηStar}
,Frame->True
,Axes->False
,FrameLabel->{"η","f(η)"}
]
]
MakeLogDensityPlot[usePars_,f_,{ηInit_,ηStar_}]:=Module[{},
LogPlot[{f[η]},{η,ηInit+SafetyEpsilon,ηStar}
,PlotRange->All
,Frame->True
,Axes->False
,FrameLabel->{"η","f(η)"}
]
]
MakeqPlot[usePars_,q_,{ηInit_,ηStar_}]:=Module[{},
Plot[{q[η]},{η,ηInit,ηStar}
,Frame->True
,Axes->False
,FrameLabel->{"η","q(η)"}
]
];
MakeiPlot[usePars_,q_,{ηInit_,ηStar_}]:=Module[{},
Plot[
Evaluate[{((-1+δ θ-θ μ+q[η])/θ),δ,δ+ρ}/.usePars]
,{η,ηInit,ηStar}
,PlotStyle->{{ColorData[1][1]},{ColorData[1][1],Dashed},{ColorData[1][1],Dashed}}
,Frame->True
,Axes->False
,FrameLabel->{"η","i(η)"}
,Epilog->{
Text["Zero growth, i == δ",{ηInit,δ/.usePars},{-1,-1}]
,Text["Growth == discount, i - δ == ρ",{ηInit,δ+ρ/.usePars},{-1,1}]
}
]
];
MakeWAPlot[usePars_,W_,WA_,{ηInit_,ηStar_}]:=Module[{},
Show[
Plot[{WA[η]},{η,ηInit,ηStar}
,PlotStyle->{{Gray}}
,Frame->True
,Axes->False
,FrameLabel->{"η","W(η)"}
],
Plot[{W[η]},{η,ηInit,ηStar}
,Frame->True
,Axes->False
,FrameLabel->{"η","W(η)"}
]
]
];
MakeWPlot[usePars_,W_,{ηInit_,ηStar_}]:=Module[{},
Plot[{W[η],W'[η]},{η,ηInit,ηStar}
,Frame->True
,Axes->False
,FrameLabel->MakeLabels["η",{"W(η)","W'(η)"}]
]
];
MakeRiskAversionPlot[usePars_,W_,q_,{ηInit_,ηStar_}]:=Module[{eps=(10^(-3))},
LogPlot[{((q'[η])/(q[η]+η))},{η,ηInit+eps,ηStar}
,Frame->True
,Axes->False
,FrameLabel->{"η","-W''(η)/W'(η)"}
]
];
MakeRecapPlot[usePars_,W_,{ηInit_,ηStar_}]:=Module[{χ=χ/.usePars,Ws=W[ηStar]},
Plot[{Ws-W[η]-(ηStar-η)-χ},{η,ηInit,ηStar}
,Frame->True
,Axes->False
,FrameLabel->{"η","W((SuperStar[η])) - W(η) - ((SuperStar[η]) - η) - χ"}
,PlotRange->All
] 
];
MakeSmallqPrimePsiPlot[usePars_,q_,{ηInit_,ηTilde_,ηStar_}]:=Module[{},
Plot[{q'[η],Clip[Ψ[η,q[η]],{0,1}]/.usePars},{η,ηInit,ηStar}
,Frame->True
,Axes->False
,FrameTicks->{{Table[i,{i,0,10}],Table[{i,""},{i,0,10}]},{{{ηInit,NumberForm[ηInit,3]},{ηTilde,NumberForm[ηTilde,3]},{ηStar,NumberForm[ηStar,3]}},None}}
,FrameLabel->MakeLabels[None,{"q'(η)","ψ(η)"}]
,Prolog->{White,Rectangle[Scaled@{0,0},Scaled@{1,1}]}
]
]
PlotTypes={"Off","Probability density","Log probability density","Pdf and cdf","Pdf and cdf w/ analytic","Pdf and cdf, no rental","Price function","Price function with inset","Price function, no rental","Value function","Value function w/ analytic","Value function, no rental","Gain from recapitalization","Gain from recapitalization, no rental","Investment rate","Investment rate, no rental","Risk aversion","Risk aversion, no rental"};
MakeAPlot[type_,usePars_,q_,f_,F_,W_,fA_,FA_,WA_,ηI_,ηT_,ηS_,qz_,fRz_,FRz_,Wz_,ηIz_,ηSz_,status_,statusz_]:=Module[{plot},
plot=Switch[type,
"Probability density",MakeDensityPlot[usePars,f,{ηI,ηS}],
"Log probability density",MakeLogDensityPlot[usePars,f,{ηI,ηS}],
"Pdf and cdf",MakePDFCDFPlot[usePars,f,F,{ηI,ηS}],
"Pdf and cdf w/ analytic",MakePDFCDFAPlot[usePars,f,F,fA,FA,{ηI,ηS}],
"Price function",MakeqPlot[usePars,q,{ηI,ηS}],
"Price function with inset",Show[MakeqPlot[usePars,q,{ηI,ηS}],Epilog->Scale[Inset[MakeSmallqPrimePsiPlot[usePars,q,{ηI,ηT,ηS}],Scaled[{0.95,0.6}],Scaled[{1,1}]],0.85]],
"Value function",MakeWPlot[usePars,W,{ηI,ηS}],
"Value function w/ analytic",MakeWAPlot[usePars,W,WA,{ηI,ηS}],
"Risk aversion",MakeRiskAversionPlot[usePars,W,q,{ηI,ηS}],
"Risk aversion, no rental",MakeRiskAversionPlot[usePars,Wz,qz,{ηI,ηS}],
"Gain from recapitalization",MakeRecapPlot[usePars,W,{ηI,ηS}],
"Gain from recapitalization, no rental",MakeRecapPlot[usePars,Wz,{ηIz,ηSz}],
"Investment rate",MakeiPlot[usePars,q,{ηI,ηS}],
"Investment rate, no rental",MakeiPlot[usePars,qz,{ηIz,ηSz}],
"Pdf and cdf, no rental",Plot[{fRz[η],FRz[η]},{η,ηIz,ηSz},Frame->True,Axes->False,FrameLabel->MakeLabels["η",{"f(η)","F(η)"}]],
(*"NetFlowReflectingBoundary",
Plot[{-ηDriftFun[η,qz[η]]fRz[η]+ηDiffusionFun[η,qz[η]]fRz'[η]}/.usePars,{η,ηIz,ηSz}
,Frame->True,Axes->False,FrameLabel->{"η","-(μ^η)f(η)+((σ^η)^2)f'(η)/2"},PlotRange->All],
"NetFlowAbsorbingBoundary",
Plot[{-ηDriftFun[η,qz[η]]fAz[η]+ηDiffusionFun[η,qz[η]]fAz'[η]}/.usePars,{η,ηIz,ηSz}
,Frame->True,Axes->False,FrameLabel->{"η","-(μ^η)f(η)+((σ^η)^2)f'(η)/2"},PlotRange->All],*)
"Price function, no rental",Plot[{qz[η]},{η,ηIz,ηSz},Frame->True,Axes->False,FrameLabel->{"η","q(η)"}],
"Value function, no rental",Plot[{Wz[η]},{η,ηIz,ηSz},Frame->True,Axes->False,FrameLabel->{"η","W(η)"}],
_,Graphics@Text["Invalid choice"]
];
plot
]
CheckConstraints[usePars_]:=Module[{errorstr="",isok=True,constraints},
constraints={
{(ρ+δ)<a<(ρ+δ+(0.5*theta*rho*rho)),"ρ < a - δ < ρ + δ + θ (ρ^2) / 2"},
{ρ>r,"ρ > r"},
{a>ab,"a > (OverBar[a])"},
{Abs[σ1]>0.0, "|σ| > 0"},
{Abs[σ1]>0.0||Abs[σ2]>0.0, "|(Subscript[σ, 1])| > 0 ∨ |(Subscript[σ, 2])| > 0"}
}/.usePars;
Do[
If[!constraint[[1]],
isok=False;
errorstr=errorstr<>"\n"<>ToString[TraditionalForm[constraint[[2]]]]
],{constraint,constraints}
];
{isok,errorstr}
];
(***** MASTER MANIPULATION FUNCTION ****)
statusThrown=-1;
MessagePanel[stuff_,title_]:=Panel[stuff,Alignment->{Top,Left},Background->Directive[Opacity[0.5],Red]]
AddToEpilog[plot_,newepi_]:=Show[plot,Epilog->Flatten[{newepi,Cases[plot,(Epilog->x_)->x,Infinity]}]];
MakeManipulatedPlots[plotTypes_,displayParameters_,ignoreConstraints_,matchqb_,usePars_]:=Module[{status=statusOK,q,f,F,fA,FA,W,WA,α,β,ηI,ηC,ηT,ηS,ηF,fAnorm,WAnorm,plots={},needRent=True,needNoRent=True,error,isok,qz,fRz,FRz,Wz,ηIz,ηSz,statusz=statusOK,WPanel,WPanelz,warns,fplotTypes,rentPlotTypes,noRentPlotTypes,make,errormsg},
If[!ignoreConstraints,
{isok,error}=CheckConstraints[usePars];
If[!isok,Return@{False,Graphics@Text["The following constraints are not satisfied: "<>error]}];
];
fplotTypes=Select[plotTypes,#≠"Off"&];
noRentPlotTypes=Select[fplotTypes,!StringFreeQ[#,"no rental"]&];
rentPlotTypes=Select[fplotTypes,StringFreeQ[#,"no rental"]&];
needNoRent=Length[noRentPlotTypes]>0;
needRent=Length[rentPlotTypes]>0;
If[needRent,
Catch[
{status,ηI,ηC,ηT,ηS,q,f,F,W}=Quiet@Check[SolveEverything[usePars],Throw[{"An unexpected message was caught -- possible reason is parameters outside the range considered typical"},"QSolver"]];
α=(1/(1+ComputeqPrimeInit[usePars]));
β=ComputefPowerLawExp[usePars];
ηF=ηI+0.001 (ηT-ηI);
fAnorm=((f[ηF])/(((ηF-ηI))^β));
fA=Function[η,fAnorm(((η-ηI))^β)];
FA=Function[η,(fAnorm/(β+1))(((η-ηI))^(β+1))];
WAnorm=((W[ηF])/(((ηF-ηI))^α));
WA=Function[η,WAnorm(((η-ηI))^α)];
,"QSolver"
,Function[{value,tag}
,status=statusThrown
;errormsg="A solution could not be reliably computed for the given parameters: "<>Fold[#1<>"\n"<>#2&,"",value]
]
];
];(* If[needRent,..] *)
If[needNoRent,
Catch[
{statusz,ηIz,ηSz,qz,fRz,FRz,Wz}=Quiet@Check[SolveEverythingNoRent[usePars],Throw[{"An unexpected message was caught -- possible reason is parameters outside the range considered typical"},"QSolver"]]
,"QSolver"
,Function[{value,tag}
,statusz=statusThrown
;errormsg="A solution could not be reliably computed for the given parameters: "<>Fold[#1<>"\n"<>#2&,"",value]
]
];
]; (* If[needNoRent,..] *)
make[type_]:=Module[{plot,stat,mods},
If[StringFreeQ[type,"no rental"]
,stat=status;mods="Model with rental"
,stat=statusz;mods="Model without rental"
];
If[stat==statusThrown
,plot=Plot[I,{x,0,1},Frame->True,Epilog->Inset[MessagePanel[errormsg,mods],Scaled@{0.5,0.5},{0,0},Scaled@{1,1}]]
,plot=MakeAPlot[type,usePars,q,f,F,W,fA,FA,WA,ηI,ηT,ηS,qz,fRz,FRz,Wz,ηIz,ηSz,status,statusz]
;If[stat≠statusOK,plot=AddToEpilog[plot,Inset[MessagePanel[StatusMsg@stat,mods],Scaled@{0.95,0.05},{1,-1},Scaled@{1,1}]]]
];
plot
];
plots=make/@fplotTypes;
If[plots=={}
,plots={If[displayParameters,Graphics@MakeParameterValues[usePars],Graphics@Text@"All plots turned off"]}
,If[displayParameters,plots[[1]]=AddToEpilog[plots[[1]],MakeParameterValues[usePars]]];
];
{status==statusOK&&statusz==statusOK,GraphicsGrid[Flatten[{{#,SpanFromLeft},{SpanFromAbove,SpanFromBoth}}&/@plots,1](*~Join~{warns}*)]}
]

plotsok=True;
{ρY,rY,σ1Y,σ2Y,aY,abY,θY,δY,χY,μY,ϵCutY,qBarSearchNDivY}={ρ,r,σ1,σ2,a,ab,θ,δ,χ,μ/δ,ϵCut,qBarSearchNDiv}/.defaultParametersRule;
Manipulate[
Module[{plots},
{plotsok,plots}=MakeManipulatedPlots[
{plot1Type,plot2Type},displayParameters,ignoreConstraints,True,
{ρ->ρX,r->rX,σ->σ1X,σ1->σ1X,σ2->σ2X,a->aX,ab->abX,θ->θX,δ->δX,χ->χX,μ->μX δX,ϵCut->ϵCutX,qBarSearchNDiv->qBarSearchNDivX}
];
If[plotsok,ρY=ρX;rY=rX;σ1Y=σ1X;σ2Y=σ2X;aY=aX;abY=abX;θY=θX;δY=δX;χY=χX;μY=μX;ϵCutY=ϵCutX;qBarSearchNDivY=qBarSearchNDivX;];
plots
]
,Style["Plot options",Bold]
,{{plot1Type,"Pdf and cdf",Tooltip["Plot 1","Functions to plot in first of two possible figures"]},PlotTypes}
,{{plot2Type,"Price function", Tooltip["Plot 2","Functions to plot in second of two possible figures"]},PlotTypes}
,{{displayParameters,True,Tooltip["Show parameters","Show selected parameter values in the figures"]},{True,False}}
,Delimiter
,Style["Model parameters",Bold]
,"Hover cursor over the parameter to get its meaning"
,{{aX,(a/.defaultParametersRule),Tooltip["a","Rate of cash production for capital at hand"]},0.0,0.3}
,{{abX,(ab/.defaultParametersRule),Tooltip["((OverBar[a])^((*)))","Rate of cash production for rented capital (only in model with option to rent)"]},0.0,0.3}
,{{δX,(δ/.defaultParametersRule),Tooltip["δ","Capital depreciation rate"]},0.0,0.3}
,{{rX,(r/.defaultParametersRule),Tooltip["r","Interest rate"]},0.01,0.1}
,{{ρX,(ρ/.defaultParametersRule),Tooltip["ρ","Discount rate"]},0.01,0.1}
,{{σ1X,(σ1/.defaultParametersRule),Tooltip["(Subscript[σ, 1])","Cash gains noise standard deviation per capital at hand"]},0.0,1.0}
,{{σ2X,(σ2/.defaultParametersRule),Tooltip["((Subscript[σ, 2])^( (*)))","Capital gains noise standard deviation per capital at hand"]},0.0,1.0}
,{{θX,(θ/.defaultParametersRule),Tooltip["θ","Curvature of the cost of adjustment"]},0.0,200.0}
,{{χX,(χ/.defaultParametersRule),Tooltip["χ","Fixed cost of recapitalization"]},0.0001,2.0}
,{{μX,(μ/δ/.defaultParametersRule),Tooltip["μ/δ","Cost of adjustment parameter"]},0.0,1.0}
(*,{{qbX,(qb/.defaultParametersRule),Tooltip["((OverBar[q])^((**)))","Initial q; -(OverBar[q]) is maximum level of debt"]},0.0,1.0}*)
(*,{{matchqb,True,Tooltip["Match initial (values^( (**) ))","Set (OverBar[q]) to the value used in the model with option to rent"]},{True,False}}*)
,Button[Tooltip["Restore parameters","If an error or warning occured, press this button to restore the most recent working parameter values"],{ρX,rX,σ1X,σ2X,aX,abX,θX,δX,χX,μX,ϵCutX,qBarSearchNDivX}={ρY,rY,σ1Y,σ2Y,aY,abY,θY,δY,χY,μY,ϵCutY,qBarSearchNDivY}
,Enabled->Dynamic@(!plotsok)
]
,Delimiter
,{{ignoreConstraints,False,Tooltip["Ignore constraints","Attempt to generate and plot solutions even if constraints are not satisfied"]},{True,False}}
,Delimiter
,"Expert parameters"
,{{ϵCutX,ϵCut/.defaultParametersRule,Tooltip["((Subscript[ϵ, xover])^( (*)))","Determines the cross-over point from approximate analytic to numerical solution near the left-hand side boundary in the model with option to rent"]}}
,{{qBarSearchNDivX,qBarSearchNDiv/.defaultParametersRule,Tooltip["(Subscript[δq, init])","Initial q search mesh subdivisions."]}}
,Delimiter
,"(*) Applies only to model with option to rent"
(*,"( * ) Applies only to model with option to rent\n ( * * ) Applies only to model without option to rent"*)
,SaveDefinitions->True
]

Publication Figures

Figures 1 and 2

This cell creates Figure 1 and Figure 2, pages 24-25 of the Nov 2012 paper

RentPlotOnePars=defPars;
RentPlotOne=SolveAndPlotEverything[RentPlotOnePars]
Export[PathToGraphics<>"RentPlotOne.pdf",RentPlotOne];
Export[PathToGraphics<>"RentPlotOne.eps",RentPlotOne];
Export[PathToGraphics<>"RentPlotOnePars.txt",RentPlotOnePars];

NoRentPlotOnePars=defPars;
NoRentPlotOne=SolveAndPlotEverythingNoRent[RentPlotOnePars]
Export[PathToGraphics<>"NoRentPlotOne.pdf",NoRentPlotOne];
Export[PathToGraphics<>"NoRentPlotOne.eps",NoRentPlotOne];
Export[PathToGraphics<>"NoRentPlotOnePars.txt",NoRentPlotOnePars];

Risk plot

This cell (a) solves everything (b) creates Figure 5 on Page 36 of the Nov 2012 version of the paper, comparing induced risk aversion -W’’/W’ for the two versions of the model

tPars=defPars;
{statusx,ηIx,ηCx,ηTx,ηSx,qnx,fnx,Fnx,Wnx}=SolveEverything[tPars];
{statusy,ηIy,ηSy,qny,fny,Fny,Wny}=SolveEverythingNoRent[tPars];
{ηIx,ηSx};
{ηIy,ηSy};

RPRent=Plot[-Wnx''[η]/Wnx'[η],{η,ηIx,ηSx}
,FrameLabel->{"η","-W''(η)/W'(η)"}
,Frame->True
,Axes->False
,Exclusions->None
,LabelStyle->12
,PlotRange->{0,4}];
RPNoRent=Plot[1/(qny[η]+η)qPrimeNoRentFun[η,qny[η]]/.tPars,{η,ηIy,ηSy}
,Frame->True
,Axes->False
,LabelStyle->12
,PlotStyle->{Dashed}
(*,PlotRange->{0.001,100}*)];
Show[RPRent,RPNoRent]
Export[PathToGraphics<>"RiskPlot.pdf",%];
Export[PathToGraphics<>"RiskPlot.eps",%%];
Export[PathToGraphics<>"RiskPlotParams.txt",tPars];

Median contour plots

These cells generate the countour plot of the median of the probability distribution. The first cell contains some adjustable parameters. MCPColorFunc determines the color function used to draw the surface; MCPNPoints determines the number of grid points to use, the higher the number, the more refined the mesh, and better the output. 

Note: A contour plot with a dense grid will take several minutes or even a better part of an hour to compute, and can consume significant amounts of memory! For quick runs, set MCPNPoints to a small value, say  ~5. Even then be prepared to wait for a few minutes. For publication quality, use MCPNPoints of about 15 - 25.

MCPColorFunc=Function[z,ColorData["SunsetColors"][1-(((1-z))^1.2)]];
MCPColorFunc=Function[z,If[z>0.0001,ColorData["RedBlueTones"][z],White]];

MCPLegend=RasterContourPlot[y,{x,0,1},{y,0,1}
,AspectRatio->8
,Contours->Table[i/10,{i,0,10}]
,ColorFunction->MCPColorFunc
,FrameTicks->{{None,Table[{i/10.0," "<>ToString@NumberForm[i/10.0,{1,1}]},{i,0,10,2}]},{None,None}}
,PlotRangePadding->{0.16,0.02}
,FrameLabel->{{None,"(OverTilde[m])"},{None,None}}
,RasterSize->{100,800}
,Background->White
];

MCPContours=Table[i/10,{i,0,10}];
MCPNPoints=10;

MedianSigmaχPlot=ContourPlot[
ComputeScaledMedian[overrideDefaultParameters[{σ1->x,χ->y}]]
,{x,0.075,0.7}
,{y,0.05,1.0}
,FrameLabel->{"σ","χ"}
,Contours->MCPContours
,ColorFunctionScaling->False
,ColorFunction->MCPColorFunc
,LabelStyle->18
,PlotRangeClipping->False
,ImageMargins->None
,ImagePadding->{{55,80},{55,10}}
,Epilog->Inset[Show[MCPLegend,LabelStyle->12],Scaled@{1.15,0.5},Scaled@{1,0.5},Scaled@{0.8,0.8}]
];
Export[PathToGraphics<>"MCPSigmaChi.pdf",MedianSigmaχPlot]

Plots for geometric proof / Ponzi illustration

RectangleIntersection[{rx0_,ry0_},{rx1_,ry1_},{px_,py_},{nx_,ny_}]:=Module[{vi,hi},
vi=Select[If[ny≠0,{#,py-(#-px)(nx/ny)}&/@{rx0,rx1},{}],ry0≤#[[2]]≤ry1&];
hi=Select[If[nx≠0,{px-(#-py)(ny/nx),#}&/@{ry0,ry1},{}],rx0<#[[1]]<rx1&];
Sort[Flatten[{vi,hi},1],#1[[1]]<#2[[1]]&]
]
Hactangle[{x0_,y0_},{x1_,y1_},α_,{d1_,d2_}]:=Module[{β,nx,ny,pts1,pts2,i,origin,oppos,left,right,line},
β=Mod[α,π];
If[β>π/2
,origin={x1,y0};oppos={x0,y1};left={x0,y0};right={x1,y1}
,origin={x0,y0};oppos={x1,y1};left={x0,y1};right={x1,y0}
];
nx=Cos[β];
ny=Sin[β];
Reap[For[i=0,True,i++,
line={};
pts1=RectangleIntersection[{x0,y0},{x1,y1},origin+i (d1+d2) {nx,ny},{nx,ny}];
pts2=RectangleIntersection[{x0,y0},{x1,y1},origin+(i (d1+d2)+d1)  {nx,ny},{nx,ny}];
If[pts1=={},Break[]];
If[pts2=={},pts2={oppos}];
line={First@pts1};
If[!(pts1[[1,1]]==pts2[[1,1]]||pts1[[1,2]]==pts2[[1,2]]),AppendTo[line, left]];
line=Join[line,pts2];
If[!(line[[-1,1]]==pts1[[-1,1]]||line[[-1,2]]==pts1[[-1,2]]),AppendTo[line, right]];
AppendTo[line,Last@pts1];
Sow[Polygon[line]];
]][[2,1]]
]

ponziPars=overrideDefaultParameters[{θ->6.5}];
noPonziPars=overrideDefaultParameters[{}]

{statusx1,ηIx1,ηCx1,ηTx1,ηSx1,qnx1,fnx1,Fnx1,Wnx1}=SolveEverything[noPonziPars];
{statusx2,ηIx2,ηCx2,ηTx2,ηSx2,qnx2,fnx2,Fnx2,Wnx2}=SolveEverything[ponziPars];
MakePonziPlot[{ηI_,ηS_,qFun_},pars_,label_:""]:=Module[{ηPPi,ηPPf,qPPi,qPPf,ηSM,qSM,noPonzi,labelSize=14},
{ηSM,qSM}=DPCTip/.pars;
ηPPi=Min[ηI,ηSM];
ηPPf=ηS;
ηPPi=ηPPi-0.1(ηPPf-ηPPi);
ηPPf=ηPPf+0.2(ηPPf-ηPPi);
qPPi=qFun[ηI]-0.1(qFun[ηS]-qFun[ηI]);
qPPf=Max[qStarsOfη/.pars/.η->ηPPf];
noPonzi=ηI>ηSM;
Show[
RegionPlot[
{Evaluate[DPCImplicitEq/.A_==0->A<0/.pars]}
,{η,ηPPi,ηPPf},{q,qPPi,qPPf}
,FrameLabel->{"η","q"}
,Prolog->{
(*{GrayLevel[0.85],Polygon[{{ηPPi,qSM},{ηPPf,qSM},{ηPPf,qPPf},{ηPPi,qPPf}},VertexTextureCoordinates->{{0,0},{1,0},{1,1},{0,1}}]}*)
{GrayLevel[0.9],Hactangle[{ηPPi,qSM},{ηPPf,qPPf},ArcTan[qPPf-qPPi,ηPPf-ηPPi],0.2(((ηPPf-ηPPi))^(-1)){0.5,0.5}]}
,{Thin,Line[{{ηPPi,qSM},{ηSM,qSM}}]}
}
,PlotRangeClipping->False
,Epilog->{
{Dashed,Gray,Line[{{ηI,qPPi},{ηI,qPPf}}]}
,{Dotted,Gray,Line[{{ηSM,qPPi},{ηSM,qPPf}}]}
,Translate[Rotate[Text[Style["η = (OverBar[η])",labelSize],{ηI,If[noPonzi,((3qPPi+qPPf)/4),((2qPPi+2qPPf)/4)]},{0,0}],90Degree],{0.03(ηPPf-ηPPi),0}]
,Translate[Rotate[Text[Style["η = (SuperStar[(Subscript[η, min])])",labelSize],{ηSM,If[noPonzi,((2qPPi+2qPPf)/4),((3qPPi+qPPf)/4)]},{0,0}],90Degree],{0.03(ηPPf-ηPPi),0}]
,Text[Style["q' < 0",labelSize],{If[noPonzi,((ηSM+2ηPPf)/3),((ηSM+ηPPf)/2)],qSM}]
,Text[Style["Growth > Discount",labelSize],If[noPonzi,{((3ηPPi+ηPPf)/4),qPPf},{((ηI+ηSM)/2),((qSM+qPPf)/2)}],{0,If[noPonzi,2,0]}]
(*,Text["Growth > Discount",{((ηI+ηSM)/2),qSM},{0,-2.25}]*)
,Point[{ηI,qFun[ηI]}]
,Point[{ηS,qFun[ηS]}]
,Text[Style[label,labelSize+2],Scaled@{-0.13,1.02},{0,1}]
}
]
,Plot[qFun[η],{η,ηI,ηS},PlotStyle->Black]
,LabelStyle->labelSize
,ImageSize->300
]
]

ponziPlotA=MakePonziPlot[{ηIx1,ηSx1,qnx1},noPonziPars,"(a)"]
ponziPlotB=MakePonziPlot[{ηIx2,ηSx2,qnx2},ponziPars,"(b)"]

Export[PathToGraphics<>"PonziPlotA.pdf",ponziPlotA]
Export[PathToGraphics<>"PonziPlotB.pdf",ponziPlotB]
Export[PathToGraphics<>"PonziPlotAParams.txt",noPonziPars]
Export[PathToGraphics<>"PonziPlotBParams.txt",ponziPars]

Probability density function plots

These cells generate the probability density plots of Fig. 3. These require some manual intervention, if the parameters are changed from baseline. Follow the instructions in the comments.

Note: These plots are somewhat slow as well, so please be patient. Rendering should not take more than a minute, though.

(* Begin by computing the critical χ values, ie. the points where recapitalization is no longer possible *)
(* Set baseline parameters but with χ now unspecified (it will be varied) *)
freeχDefPars=overrideDefaultParameters[{χ->χ}]
(* Since we only need a plot or two of this kind, it is not worth the effort to make everything fully automatic. Some manual intervention is therefore needed, namely, one should set PDPχBegin and PDPχEnd parameters so that that range brackets a root of recapitalization gain for the chosen set of parameters (experiment with the plots below; currently works for default parameters). Below one also needs to enter a guess for the critical χ values manually. *)
PDPχBegin=0.5; (* <-- Serves as the lower bound on χ in the following plots *)
PDPχEnd=0.9; (* <-- This is the corresponding upper bound *)
PDPnχ=200; (* <-- Sets the number of χ values on which the problem is actually solved. Set this to low values (~10) for quick and dirty plots; Something ~100 is recommended for production quality figures. *)
PDPnη=100;  (* <-- Sets the number of η values for which the density is evaluated *)
Plot[ComputeRecapitalizationGainNoRent[-ηSingular/.freeχDefPars,freeχDefPars],{χ,PDPχBegin,PDPχEnd}]
Plot[ComputeRecapitalizationGain[-ηSingular/.freeχDefPars,freeχDefPars],{χ,PDPχBegin,PDPχEnd}]

(* The above two plots should be now inspected, and a rough estimate for a root is then entered as values of PDPχGuessNR and PDPχGuess, corresponding to models without and with rent, respectively. One should make sure that the output of the FindRoots on the two last lines is sane. *)
PDPχGuessNR=0.75;
PDPχGuess=0.75;
χCritNR=χ/.FindRoot[ComputeRecapitalizationGainNoRent[-ηSingular/.freeχDefPars,freeχDefPars][[2]],{χ,PDPχGuessNR}]
χCrit=χ/.FindRoot[ComputeRecapitalizationGain[-ηSingular/.freeχDefPars,freeχDefPars][[2]],{χ,PDPχGuess}]

(* Let us next generate the values of the distribution function on a regular grid. We could simply construct a function that solves the problem for any input (χ, η), but that would lead to mathematica solving the problem completely from scratch for each value of η -- a massive waste of cpu cycles, and time! So, let's just tabulate the densities. The mesh density (and thus the level of detail in the plot) is controlled by PDPnχ and PDPnη parameters, set earlier. The higher these values, the nicer the plot (and slower the rendering). *)
PDPηMin=100.0; (* <-- min. value of η seen in data, init'd to an arb. large value *)
PDPηMax=-100.0; (* <-- max. value of η seen in data, init'd to an arb. small value *)
PDPDataNR=Reap[
Do[
{statusx,ηIx,ηSx,qnx,fnx,Fnx,Wnx}=SolveEverythingNoRent[freeχDefPars];
PDPηMin=Min[ηIx,PDPηMin];
PDPηMax=Max[ηSx,PDPηMax];
Do[With[{η=ηIx+If[iη==0,1.0 (10^(-3))(ηSx-ηIx),((ηSx-ηIx)/(PDPnη))iη]},Sow@{η,χ,fnx[η]}],{iη,0,PDPnη}
],{χ,PDPχBegin,PDPχEnd,((PDPχEnd-PDPχBegin)/(PDPnχ))}
]
][[2,1]];
PDPData=Reap[
Do[
{statusx,ηIx,ηCx,ηTx,ηSx,qnx,fnx,Fnx,Wnx}=SolveEverything[freeχDefPars];
PDPηMin=Min[ηIx,PDPηMin];
PDPηMax=Max[ηSx,PDPηMax];
Do[With[{η=ηIx+If[iη==0,1.0 (10^(-3))(ηSx-ηIx),((ηSx-ηIx)/(PDPnη))iη]},Sow@{η,χ,fnx[η]}],{iη,0,PDPnη}
],{χ,PDPχBegin,PDPχEnd,((PDPχEnd-PDPχBegin)/(PDPnχ))}
]
][[2,1]];
(* Specially solve the distribution on the critical χ value -- just to draw a curve on the surface to indicate the cross-over χ from recapitalization to liquidation *)
{PDPCritηstatusx,PDPCritηIx,PDPCritηCx,PDPCritηTx,PDPCritηSx,PDPCritqnx,PDPCritfnx,PDPCritFnx,PDPCritWnx}=SolveEverything[overrideDefaultParameters[{χ->χCrit}]];
{PDPCritNRηstatusx,PDPCritNRηIx,PDPCritNRηSx,PDPCritNRqnx,PDPCritNRfnx,PDPCritNRFnx,PDPCritNRWnx}=SolveEverythingNoRent[overrideDefaultParameters[{χ->χCritNR}]];

(* Function to create custom ticks and tick labels *)
PDPTicks[xlo_,xhi_,unit_,subdiv_,fun_]:=Module[
{i0=Floor[subdiv xlo/unit],
i1=Ceiling[subdiv xhi/unit]},
Table[If[Mod[i,subdiv]==0,{(i/subdiv)unit,fun[(i/subdiv)unit],{0.01,0}},{(i/subdiv)unit,"",{0.005,0}}],{i,i0,i1}]
];
(* Main plot in the case with rental *)
PDPlot=Show[
ListPlot3D[PDPData
,PlotRange->{Automatic,Automatic,{0,3}}
,BoxRatios->{1,1,1}
,AxesLabel->{"η","      χ","f"}
,LabelStyle->16
,Mesh->10
,Ticks->{
PDPTicks[-1,1.5,0.5,5,NumberForm[#,{2,1}]&],
PDPTicks[0.2,0.8,0.1,5,NumberForm[#,{2,1}]&],
PDPTicks[0,3,1,5,NumberForm[#,1]&]
}
,BoxStyle->Directive[Black,Thickness[0.0005]]
,ViewPoint->{5,-10,7.5}
,AxesStyle->Thickness[0.0025]
],
ParametricPlot3D[{η,χCrit,PDPCritfnx[η]},{η,PDPCritηIx,PDPCritηSx},PlotPoints->30,PlotStyle->{{Black,Thick}}],
Graphics3D@{Black,Dashed,Line[{{PDPηMin,χCrit,0.0},{PDPηMax,χCrit,0.0}}]},
Graphics3D@{Text[Style["(OverBar[χ])",14],{PDPηMax,χCrit,0},{-3.0,0}]}
]
(* Same, but without rental *)
PDPlotNR=Show[
ListPlot3D[PDPDataNR
,PlotRange->{Automatic,Automatic,{0,3}}
,BoxRatios->{1,1,1}
,AxesLabel->{"η","      χ","f"}
,LabelStyle->16
,Mesh->10
,Ticks->{
PDPTicks[-1,1.5,0.5,5,NumberForm[#,{2,1}]&],
PDPTicks[0.2,0.8,0.1,5,NumberForm[#,{2,1}]&],
PDPTicks[0,3,1,5,NumberForm[#,1]&]
}
,BoxStyle->Directive[Black,Thickness[0.0005]]
,ViewPoint->{5,-10,7.5}
,AxesStyle->Thickness[0.0025]
],
ParametricPlot3D[{η,χCritNR,PDPCritNRfnx[η]+0.01},{η,PDPCritNRηIx,PDPCritNRηSx},PlotPoints->30,PlotStyle->{{Black,Thick}}],
ParametricPlot3D[{η,χCritNR,PDPCritNRfnx[η]-0.01},{η,PDPCritNRηIx,PDPCritNRηSx},PlotPoints->30,PlotStyle->{{Black,Thick}}],
Graphics3D@{Black,Dashed,Line[{{PDPηMin,χCritNR,0.0},{PDPηMax,χCritNR,0.0}}]},
Graphics3D@{Text[Style[" (OverBar[χ])",14],{PDPηMax,χCritNR,0},{-5.0,0}]}
]
Export[PathToGraphics<>"PDPlot.pdf",Rasterize[PDPlot,RasterSize->4000]]
Export[PathToGraphics<>"PDPlot.eps",Rasterize[PDPlot,RasterSize->4000]]
Export[PathToGraphics<>"PDPlotNR.pdf",Rasterize[PDPlotNR,RasterSize->4000]]
Export[PathToGraphics<>"PDPlotNR.eps",Rasterize[PDPlotNR,RasterSize->4000]]
Export[PathToGraphics<>"PDPlotParameters.txt",freeχDefPars]