Post-crisis macrofinancial modelling: continuous time approaches

Isohätälä, Klimenko, and Milne - Illustrative example numerical solver

This is a notebook for the numerical solution and presentation of some of the steps in the derivation of a model described in “The Handbook of Post Crisis Financial Modelling”, Palgrave MacMillan Handbook Series, eds. E. Haven et al., chapter “Post-crisis macrofinancial modelling: continuous time approaches”.

If interested only in the numerical output of the model, evaluate the notebook evaluation cells, and scroll down to the interactive solver at the end. The notebook initial sections contain some of the explicit steps, with commentary, needed in constructing the final equations.
Note: The numerical solver is tested and known to work well at least for parameters close to baseline. Some extreme parameter values may cause the solver to fail, or produce unreliable results without warning. We have tried our best to make the routines robust, but cannot guarantee that the output is correct for all input parameters. Do not hesitate to contact the authors for assistance.

Derivations

Here we give the model equations, beginning right from the stochastic differential equations. We state all the equations necessary to complete the model, and use Mathematica to verify that what we have written indeed is correct. These derivation steps are not necessary for the numerical solution of the model, but are useful in ensuring the correctness of our formulas.

Ito differentiation rule

The Ito's Lemma is used frequently in the derivations, and for this, we need to write our own function.

ItoD::usage="Apply the Ito's Lemma on a given expression. ItoD[F, vars, SDES, dt, dz] differentiates F with respect to the variables vars with differentials SDEs where dt is the time differential and dz are differentials of all the independent Weiner process in SDEs.";
ItoD[F_,{vars__},{SDEs__},dt_,{dz__}]:=Module[
{μ=Coefficient[{SDEs},dt],
σ=Normal@CoefficientArrays[{SDEs},{dz}][[2]]},
(((∂_({{vars}}))F).μ+(1/2)Tr[(σ^\[Transpose]).((∂_({{vars}},{{vars}}))F).σ])dt+(((∂_({{vars}}))F).σ).{dz}]

Stochastic differential equations

These form the basis of the model: We give the differentials  dc and dk and postulate some yet to be determined price and aggregate state processes q and x. From c,k variables we then move to n and ϕ, where ϕ is a choice and so it doesn’t need its own differential equation.

Dc=(a k+r c-q τ k-κ k)\[DifferentialD]t;
Dk=τ k \[DifferentialD]t+σ k \[DifferentialD]z;
Dq=μq q \[DifferentialD]t+σq q \[DifferentialD]z;
Dx=μx x \[DifferentialD]t+σx x \[DifferentialD]z;
Dn=(r+((a/q)+μq+σ σq-r)ϕ-λ) n\[DifferentialD]t+(σ+σq) ϕ n\[DifferentialD]z;

Verify formulas derived by hand: Ito differentiate definition of n and compare to what we defined above.

0== Dn-ItoD[c+q k,{c,k,q},{Dc/.κ->(λ/k)n,Dk,Dq},\[DifferentialD]t,{\[DifferentialD]z}]/.{c->(1-ϕ)n,k->ϕ n/q}//Simplify

Hamilton-Jacobi-Bellman and first order conditions

Next, we state the Hamilton-Jacobi-Bellman (HJB) equations, in full form for the value function V(n,x) and in the reduced form for the marginal value W(x), where V(n,x)=n W(x). We also state optimal choice conditions for asset allocations ϕ and (OverBar[ϕ]) and consumption λ (these are the expressions ϕZero, ϕZeroB, and λZero).

ϕZero=(a/(q[x]))+μq[x]+σ σq[x]-r+x (σ+σq[x]) σx[x] ((W'[x])/(W[x]));
λZero=1-W[x];
HJBFull=-ρ V[n,x]+λ n+(r+((a/(q[x]))+μq[x]+σ σq[x]-r)ϕ-λ)n (∂_n)V[n,x]+μx[x]x(∂_x)V[n,x]+(1/2)(((σx[x]x))^2)(∂_(x,x))V[n,x]+σx[x]x(σ+σq[x])ϕ n (∂_(n,x))V[n,x]+(1/2)(((σ+σq[x]))^2)(ϕ^2)(n^2)(∂_(n,n))V[n,x];
HJB=λ λZero+ϕ ϕZero W[x]+(r-ρ) W[x]+x μx[x] W'[x]+(1/2)(x^2) ((σx[x])^2) W''[x];
HJBZero=HJB/.{λ->0,ϕ->0};
{ϕZeroB,DnB}={ϕZero,Dn}/.{a->ab,W->Wb,ρ->r,ϕ->ϕb,λ->λb}/.Wb->(1&);

Rule to map (μ^q), (σ^q), q, (μ^x), and (σ^x) to functions of the aggregate state x.

xifyRule={f:(q|μq|σq|μx|σx):>f[x]};

Verify the above HJB equations:

HJBAuto=-ρ V[n,x]+λ n+Coefficient[ItoD[V[n,x],{n,x},{Dn,Dx},\[DifferentialD]t,{\[DifferentialD]z}],\[DifferentialD]t];
0==HJBFull-(HJBAuto/.xifyRule)//Simplify
0==HJB-((1/n)HJBAuto/.V->Function[{n,x},n W[x]]/.xifyRule)//Simplify

Aggregate state equations of motion / closure conditions

Here we derive the expressions for the q drift and diffusion as given by the aggregate state equations of motion, and state and verify the aggregate drift and diffusion.

DqFormal=ItoD[q[x],{x},{Dx/.xifyRule},\[DifferentialD]t,{\[DifferentialD]z}]-(Dq/.xifyRule);
qddRule=First@Solve[#==0&/@Thread@Coefficient[DqFormal,{\[DifferentialD]t,\[DifferentialD]z}],{μq[x],σq[x]}]//Simplify;

DKtot=σ Ktot \[DifferentialD]z;
DxReal= ((a/q)ψ+(μq-(σ^2)-σ σq-((σq)^2)-r) (ψ-x)-λ x) \[DifferentialD]t+(σ+σq) (ψ-x)\[DifferentialD]z;
ClosureZero=Thread@Coefficient[DxReal-Dx/.xifyRule,{\[DifferentialD]t,\[DifferentialD]z}]/.λ->0/.qddRule;

Verify formulas derived by hand. The above aggregate x differential should be the Ito differential of the expert net worth share.

0==DxReal-ItoD[(n/(q Ktot)),{n,q,Ktot},{Dn,Dq,DKtot},\[DifferentialD]t,{\[DifferentialD]z}]/.Ktot->(n/(q x))/.ϕ->ψ/x//FullSimplify

Closed form solution of aggregate dynamics and capital allocations

The final model differential equations, along with the aggregate state drift and diffusion, are obtained by solving the appropriate equations, namely, the optimal asset allocation conditions, closure conditions, and the expert HJB. In the definitions below, we use the short-hand definitions p(x)== x q'(x)/q(x) and w(x)=x W'(x)/W(x).

Parameters are given in vector form. A symbolic placeholder definition, with helper functions to access the parameter vector:

par={a,ab,σ,r,ρ};
pa=#[[1]]&;pab=#[[2]]&;pσ=#[[3]]&;pr=#[[4]]&;pρ=#[[5]]&;

The HJB equation has the same form regardless of whether experts hold capital or not.

WPrimePrimeFun[x_,q_,p_,W_,w_,μx_,σx_,{a_,ab_,σ_,r_,ρ_}]=-((2W)/((x^2) ((σx)^2))) (r-ρ+μx w);

In the regime (call that regime I) where households hold capital, the solution of the expert and households optimal asset allocation equation (here ϕZero = 0, ϕZeroB = 0), and the closure conditions (ClosureZero =  0), and the expert HJB (HJBZero = 0) gives the following for (σ^x), (μ^x), ψ,  q'' (W'' is given above as its form is the same whether households hold capital or not):

σxFunI[x_,q_,p_,w_,{a_,ab_,σ_,r_,ρ_}]=-((a-ab)/(q w ((σ/2)+(Sqrt[((((σ/2)))^2)-p((a-ab)/(q w))]))));
μxFunI[x_,q_,p_,w_,σx_,{a_,ab_,σ_,r_,ρ_}]=(a/q)-σx(σ+σx(w+p));
ψFunI[x_,q_,p_,w_,σx_,{a_,ab_,σ_,r_,ρ_}]=x(1+((σx)/(σ+σx p)));
qPrimePrimeFunI[x_,q_,p_,w_,σx_,{a_,ab_,σ_,r_,ρ_}]=-((2 q)/((x^2)((σx)^2))) ((a/q)(1+p)-r-σx(σx (p^2)-σ w));

In regime II, that is, when households do not hold capital, the solution of the asset allocation condition (ϕZero = 0, note that ϕZeroB < 0 so that ψ = 1),  closure conditions (ClosureZero = 0), and the expert HJB gives:

σxFunII[x_,q_,p_,w_,{a_,ab_,σ_,r_,ρ_}]=(((1-x) σ)/(x-(1-x)p));
μxFunII[x_,q_,p_,w_,σx_,{a_,ab_,σ_,r_,ρ_}]=(a/q)-((1-x)/x)(σ+σx p) (σ+σx (p+w));
qPrimePrimeFunII[x_,q_,p_,w_,σx_,{a_,ab_,σ_,r_,ρ_}]=(2/((x^2)((σx)^2)))q(r-(1+p)(a/q)-(σ+σx p)(σx w-((1-x)/x)p (σ+σx (p+w)))-σ σx p);

All the above formulas are parametrized for given (σ^x) or (μ^x), (σ^x). For conveniance, we define alternative forms with the aggregate drift and diffusion replaced by their appropriate values.

ψXFunI[x_,q_,p_,w_,{a_,ab_,σ_,r_,ρ_}]=ψFunI[x,q,p,w,σxFunI[x,q,p,w,{a,ab,σ,r,ρ}],{a,ab,σ,r,ρ}];
μxXFunI[x_,q_,p_,w_,{a_,ab_,σ_,r_,ρ_}]=μxFunI[x,q,p,w,σxFunI[x,q,p,w,{a,ab,σ,r,ρ}],{a,ab,σ,r,ρ}];
μxXFunII[x_,q_,p_,w_,{a_,ab_,σ_,r_,ρ_}]=μxFunII[x,q,p,w,σxFunII[x,q,p,w,{a,ab,σ,r,ρ}],{a,ab,σ,r,ρ}];
qPrimePrimeXFunI[x_,q_,p_,w_,{a_,ab_,σ_,r_,ρ_}]=qPrimePrimeFunI[x,q,p,w,σxFunI[x,q,p,w,{a,ab,σ,r,ρ}],{a,ab,σ,r,ρ}];
qPrimePrimeXFunII[x_,q_,p_,w_,{a_,ab_,σ_,r_,ρ_}]=qPrimePrimeFunII[x,q,p,w,σxFunII[x,q,p,w,{a,ab,σ,r,ρ}],{a,ab,σ,r,ρ}];
WPrimePrimeXFunI[x_,q_,p_,W_,w_,{a_,ab_,σ_,r_,ρ_}]=Module[{μx,σx},
σx=σxFunI[x,q,p,w,{a,ab,σ,r,ρ}];
μx=μxFunI[x,q,p,w,σx,{a,ab,σ,r,ρ}];
WPrimePrimeFun[x,q,p,W,w,μx,σx,{a,ab,σ,r,ρ}]
];
WPrimePrimeXFunII[x_,q_,p_,W_,w_,{a_,ab_,σ_,r_,ρ_}]=Module[{μx,σx},
σx=σxFunII[x,q,p,w,{a,ab,σ,r,ρ}];
μx=μxFunII[x,q,p,w,σx,{a,ab,σ,r,ρ}];
WPrimePrimeFun[x,q,p,W,w,μx,σx,{a,ab,σ,r,ρ}]
];

Finally, these are the exponents used in the asymptotic expansions of q, W, and f: q(x) = (Subscript[q, 1])(x^α)+ o((x^α)), W(x) = (Subscript[W, 1])(x^(-β))+o((x^(-β))), f(x)=(Subscript[f, 1])(x^γ)+o((x^γ)). Constants (Subscript[q, 1]), (Subscript[W, 1]), (Subscript[f, 1]) are boundary condition / normalisation dependent.

βFun[{a_,ab_,σ_,r_,ρ_}]=(1/2) (1-((ab ρ)/(a r))-(((((a-ab))^2) r)/(2 a ab (σ^2)))+(Sqrt[(((1-((ab ρ)/(a r))-(((((a-ab))^2) r)/(2 a ab (σ^2)))))^2)+((2 (((a-ab))^2) r)/(a ab (σ^2)))]));
αFun[{a_,ab_,σ_,r_,ρ_}]=(1/2)(1-2 β (1+((a ab β (σ^2))/((((a-ab))^2) r)))+(Sqrt[(((1-2 β (1+((a ab β (σ^2))/((((a-ab))^2) r)))))^2)+((8 ((ab^2) (β^2) (σ^2)))/((((a-ab))^2) r))]))/.β->βFun[{a,ab,σ,r,ρ}];
γFun[{a_,ab_,σ_,r_,ρ_}]=2 (-1+β+((ab (ab+a (-1+β)) β (σ^2))/((((a-ab))^2) r)))/.β->βFun[{a,ab,σ,r,ρ}];

Finally, and most importantly, we verify that our solution indeed satisfies the appropriate equations:

{0,0,0,0,0}=={HJBZero,ϕZeroB,ϕZero}~Join~ClosureZero/.qddRule//.{
W''[x]->WPrimePrimeFun[x,q[x],((x q'[x])/(q[x])),W[x],((x W'[x])/(W[x])),μx[x],σx[x],par],
q''[x]->qPrimePrimeFunI[x,q[x],((x q'[x])/(q[x])),((x W'[x])/(W[x])),σx[x],par],
σx[x]->σxFunI[x,q[x],((x q'[x])/(q[x])),((x W'[x])/(W[x])),par],
μx[x]->μxFunI[x,q[x],((x q'[x])/(q[x])),((x W'[x])/(W[x])),σx[x],par],
ψ->ψFunI[x,q[x],((x q'[x])/(q[x])),((x W'[x])/(W[x])),σx[x],par]
}//Simplify

{0,0,0,0}=={HJBZero,ϕZero}~Join~ClosureZero/.qddRule//.{
W''[x]->WPrimePrimeFun[x,q[x],((x q'[x])/(q[x])),W[x],((x W'[x])/(W[x])),μx[x],σx[x],par],
q''[x]->qPrimePrimeFunII[x,q[x],((x q'[x])/(q[x])),((x W'[x])/(W[x])),σx[x],par],
σx[x]->σxFunII[x,q[x],((x q'[x])/(q[x])),((x W'[x])/(W[x])),par],
μx[x]->μxFunII[x,q[x],((x q'[x])/(q[x])),((x W'[x])/(W[x])),σx[x],par],
ψ->1
}//Simplify

Numerical solver

Here we give the actual numerical solver. SolveQW is the routine that does the bulk of the work: For a given parameters and (Subscript[q, 1]), that is, the coefficient on the asymptotic expansion of q, it solves the q and W differential equations upto the point (SuperStar[x]) where W'((SuperStar[x]))=1. This is essentially done in three parts: For x∈[0,ϵ], the asymptotic expansion is used (ϵ is chosen as the square root of the machine epsilon). From x = ϵ upto the point where ψ=1, we use the standard Mathematica ordinary differential equation solver to solve the regime I equations, setup so that it automatically stops at first x such that ψ=1 (call that (OverTilde[x])). Then, from (OverTilde[x]) onwards to (SuperStar[x]) we do the same for regime II equations, now stopping at first x such that W'(x) = 0.

In order to fix (Subscript[q, 1]) so that q'((SuperStar[x]))=0, we define the helper function TestFun (in the text we called this Θ) that simply uses SolveQW to find q and W, and then outputs q'((SuperStar[x])). We can then just use the Mathematica standard root solver to find the corrext (Subscript[q, 1]). SolveQWAuto function attempts to do this root finding automatically. It assumes a certain search range for the correct (Subscript[q, 1]) which might not be always appropriate. If it fails, then one should use TestFun and SolveQW “manually” to pin down the right solution.

We then give the ergodic distribution f solver SolveErgodic. All this does is it solves the Kolmogorov Forward Equation / Fokker-Planck equation with given input q and W. The numerical solver does not directly solve f, but g(x)=(((σ^x)(x))^2)f(x)/2 instead. The full x range is again spliced into three intervals. Asymptotic expansion is used in the first, then Region I equations, and then Region II.

Finally, there is the auxilliary function SolveOthers, which gives the capital allocation function ψ, and the price and aggregate state drift and diffusion functions.

SolveQW[q1_?NumericQ,par:{a_,ab_,σ_,r_,ρ_}]:=Module[
{ϵ=($MachineEpsilon^(1/2))
,success=False,xtilde,xstar
,α,β
,q0n,qp0n,W0n,Wp0n
,qIn,qpIn,WIn,WpIn
,qIIn,qpIIn,WIIn,WpIIn
,qn,qpn,Wn,Wpn
,Wnorm
,qPrimeMax=1.0 (10^2)
},

(* Analytic solution from asympt. expansion *)
α=αFun[par];
β=βFun[par];
q0n=x|->(ab/r)+q1 (x^α);
qp0n=x|->α q1 (x^(α-1));
W0n=x|->(x^(-β));
Wp0n=x|->-β (x^(-β-1));

(* Region I solution, taking off from the asymptotic form at some small x = ϵ. *)
{qIn,qpIn,WIn,WpIn}={q,q',W,W'}/.Quiet@First@NDSolve[{
q''[x]==qPrimePrimeXFunI[x,q[x],((x q'[x])/(q[x])),((x W'[x])/(W[x])),par],
W''[x]==WPrimePrimeXFunI[x,q[x],((x q'[x])/(q[x])),W[x],((x W'[x])/(W[x])),par],
q[ϵ]==q0n[ϵ],q'[ϵ]==qp0n[ϵ],W[ϵ]==W0n[ϵ],W'[ϵ]==Wp0n[ϵ]}
,{q,q',W,W'}
,{x,ϵ,1}
,Method->{"EventLocator"
,"Event"->{ψXFunI[x,q[x],((x q'[x])/(q[x])),((x W'[x])/(W[x])),par]-1,q[x]-(a/r)}
,"EventAction":>{Throw[success=True;xtilde=x,"StopIntegration"],Throw[success=False,"StopIntegration"]}
}];

If[!success,Return[{success,{qIn,qpIn,WIn,WpIn},{ϵ,qIn[[1,1,2]],qIn[[1,1,2]]}}]];

(* Region II solution, continuing from where Region I solver stopped. *)
success=False;
{qIIn,qpIIn,WIIn,WpIIn}={q,q',W,W'}/.Quiet@First@NDSolve[{
q''[x]==qPrimePrimeXFunII[x,q[x],((x q'[x])/(q[x])),((x W'[x])/(W[x])),par],
W''[x]==WPrimePrimeXFunII[x,q[x],((x q'[x])/(q[x])),W[x],((x W'[x])/(W[x])),par],
q[xtilde]==qIn[xtilde],q'[xtilde]==qpIn[xtilde],
W[xtilde]==WIn[xtilde],W'[xtilde]==WpIn[xtilde]}
,{q,q',W,W'}
,{x,xtilde,1}
,Method->{"EventLocator"
,"Event"->{W'[x],q[x]-(a/r)}
,"EventAction":>{Throw[success=True;xstar=x,"StopIntegration"],Throw[success=False,"StopIntegration"]}
}];

(* Define the whole solution piecewise of the above solutions *)
Wnorm=((WIIn[WIIn[[1,1,2]]])^(-1));
qn=x|->({{\[Piecewise], ({{q0n[x], x<ϵ}, {qIn[x], ϵ≤x<xtilde}, {qIIn[x], x≥xtilde}})}});qpn=x|->({{\[Piecewise], ({{qp0n[x], x<ϵ}, {qpIn[x], ϵ≤x<xtilde}, {qpIIn[x], x≥xtilde}})}});
Wn=x|->({{\[Piecewise], ({{Wnorm W0n[x], x<ϵ}, {Wnorm WIn[x], ϵ≤x<xtilde}, {Wnorm WIIn[x], x≥xtilde}})}});Wpn=x|->({{\[Piecewise], ({{Wnorm Wp0n[x], x<ϵ}, {Wnorm WpIn[x], ϵ≤x<xtilde}, {Wnorm WpIIn[x], x≥xtilde}})}});

If[!success,Return[{success,{qn,qpn,Wn,Wpn},{0,xtilde,qIIn[[1,1,2]]}}]];

(* Successful exit *)
{success,{qn,qpn,Wn,Wpn},{0,xtilde,xstar}}
]
TestFun::usage="Solve the model for the q and W functions, and simply return the value of q' at the final point. When this is zero, all the boundary conditions are satisfied, and we have guessed the correct value of q1.";
TestFun[q1_?NumericQ,par:{a_,ab_,σ_,r_,ρ_}]:=Module[{success,q,qp,W,Wp,xzero,xtilde,xstar},
{success,{q,qp,W,Wp},{xzero,xtilde,xstar}}=SolveQW[q1,par];
qp[xstar]
]

SolveErgodic[{q_,qp_,W_,Wp_},{xtilde_,xstar_},par:{a_,ab_,σ_,r_,ρ_}]:=Module[{gI,FI,gII,FII,sI,sII,RI,RII,Nn,Foff,f0,F0,γ,ϵ=($MachineEpsilon^(1/2))},
γ=γFun[par];
(* If γ < -1, all probability mass drift towards zero. *)
If[γ<-1,
Return@{x|->DiracDelta[x],x|->UnitStep[x],{0,xstar}}
];

sI[x_?NumericQ]=(x^2)((σxFunI[x,q[x],((x qp[x])/(q[x])),((x Wp[x])/(W[x])),par])^2);
sII[x_?NumericQ]=(x^2)((σxFunII[x,q[x],((x qp[x])/(q[x])),((x Wp[x])/(W[x])),par])^2);
RI[x_?NumericQ]=2μxXFunI[x,q[x],((x qp[x])/(q[x])),((x Wp[x])/(W[x])),par]/(x ((σxFunI[x,q[x],((x qp[x])/(q[x])),((x Wp[x])/(W[x])),par])^2));
RII[x_?NumericQ]=2μxXFunII[x,q[x],((x qp[x])/(q[x])),((x Wp[x])/(W[x])),par]/(x ((σxFunII[x,q[x],((x qp[x])/(q[x])),((x Wp[x])/(W[x])),par])^2));

f0=x|->(x^γ);
F0=x|->(1/(γ+1))(x^(γ+1));
{gI,FI}={g,F}/.First@NDSolve[{0==RI[x]g[x]-g'[x],F'[x]==((2 g[x])/(sI[x])),g[ϵ]==(1/2)sI[ϵ]f0[ϵ],F[ϵ]==F0[ϵ]},{g,F},{x,ϵ,xtilde}];
{gII,FII}={g,F}/.First@NDSolve[{0==RII[x]g[x]-g'[x],F'[x]==((2 g[x])/(sII[x])),g[xtilde]==gI[xtilde],F[xtilde]==FI[xtilde]},{g,F},{x,xtilde,xstar}];

Nn=((FII[xstar])^(-1));

{x|->({{\[Piecewise], ({{Nn f0[x], x<ϵ}, {Nn((2gI[x])/(sI[x])), ϵ≤x<xtilde}, {Nn ((2gII[x])/(sII[x])), x≥xtilde}})}}),x|->({{\[Piecewise], ({{Nn F0[x], x<ϵ}, {Nn FI[x], ϵ≤x<xtilde}, {Nn FII[x], x≥xtilde}})}}),
{0,gII[[1,1,2]]}}
]

SolveQWAuto[par:{a_,ab_,σ_,r_,ρ_}]:=Module[{q0=1.0 (10^(-4)),q1=0.5,q1Max=3.0,t0,t1,step=0.5,qstar,loop=True},

t0=TestFun[q0,par];
t1=TestFun[q1,par];
While[loop,
If[t0 t1≤0.0
,qstar=qq/.FindRoot[TestFun[qq,par],{qq,0.5(q0+q1),q0,q1}]
;loop=False
,q0=q1;t0=t1
;q1=q1+step;t1=TestFun[q1,par]
];
If[q1>q1Max,Throw["Failed bracketing q1 root","IKMFail"]];
];

SolveQW[qstar,par]
]

SolveOthers[{q_,qp_,W_,Wp_},{xtilde_,xstar_},par:{a_,ab_,σ_,r_,ρ_}]:=Module[{σx,μx,ψ,qpp,μq,σq},
ψ=x|->({{\[Piecewise], ({{ψXFunI[x,q[x],((x qp[x])/(q[x])),((x Wp[x])/(W[x])),par], x<xtilde}, {1, x≥xtilde}})}});
μx=x|->({{\[Piecewise], ({{μxXFunI[x,q[x],((x qp[x])/(q[x])),((x Wp[x])/(W[x])),par], x<xtilde}, {μxXFunII[x,q[x],((x qp[x])/(q[x])),((x Wp[x])/(W[x])),par], x≥xtilde}})}});
σx=x|->({{\[Piecewise], ({{σxFunI[x,q[x],((x qp[x])/(q[x])),((x Wp[x])/(W[x])),par], x<xtilde}, {σxFunII[x,q[x],((x qp[x])/(q[x])),((x Wp[x])/(W[x])),par], x≥xtilde}})}});
qpp=x|->({{\[Piecewise], ({{qPrimePrimeXFunI[x,q[x],((x qp[x])/(q[x])),((x Wp[x])/(W[x])),par], x<xtilde}, {qPrimePrimeXFunI[x,q[x],((x qp[x])/(q[x])),((x Wp[x])/(W[x])),par], x≥xtilde}})}});
μq=x|->μx[x]((x qp[x])/(q[x]))+(1/2)((σx[x])^2)(((x^2) qpp[x])/(q[x]));
σq=x|->σx[x]((x qp[x])/(q[x]));
{qpp,ψ,μx,σx,μq,σq}
]

Auxilliary functions

These only affect the styling of the output figures.

Plot styling

labelSize=14;
smallBoost=4;

xlabel=Style["x",Italic];
Wlabel=Style[Row[{Style["W",Italic],"(",xlabel,")"}]];
Wplabel=Style[Row[{Style["W",Italic],"'(",xlabel,")"}]];
qlabel=Style[Row[{Style["q",Italic],"(",xlabel,")"}]];
qplabel=Style[Row[{Style["q",Italic],"'(",xlabel,")"}]];
ψlabel=Style[Row[{"ψ(",xlabel,")"}]];
σxlabel=Style[Row[{Superscript["σ",xlabel],"(",xlabel,")"}]];
μxlabel=Style[Row[{Superscript["μ",xlabel],"(",xlabel,")"}]];
σqlabel=Style[Row[{Superscript["σ",Style["q",Italic]],"(",xlabel,")"}]];
μqlabel=Style[Row[{Superscript["μ",Style["q",Italic]],"(",xlabel,")"}]];
flabel=Style[Row[{Style["f",Italic],"(",xlabel,")"}]];
Flabel=Style[Row[{Style["F",Italic],"(",xlabel,")"}]];
fFlabel=Style[Row[{Style["f",Italic],"(",xlabel,"),   ",Style["F",Italic],"(",xlabel,")"}]];

Plot functions

Clear[PlotBasic,PlotBasicFill,PlotBasicChoice]

Options[PlotBasic]={ShowGuides->False};
PlotBasic[{fun_,label_},{xzero_,xtilde_,xstar_},OptionsPattern[]]:=Plot[fun[x],{x,xzero,xstar}
,Axes->False,Frame->True
,FrameLabel->{xlabel,label}
,LabelStyle->labelSize
,Exclusions->{xtilde}
,AspectRatio->1
,Epilog->If[OptionValue[ShowGuides],{
{Gray,Dashed,Line[{{xtilde,-(10^4)},{xtilde,+(10^4)}}]}
,{Text["ψ < 1",Scaled@{(xtilde/xstar),0.5},{0,-2},{0,1}]}
,{Text["ψ = 1",Scaled@{(xtilde/xstar),0.5},{0,2},{0,1}]}
},{}]
,ImageSize->300
]
Options[PlotBasicFill]={ShowGuides->False};
PlotBasicFill[{fun_,label_},{xzero_,xtilde_,xstar_},OptionsPattern[]]:=Plot[fun[x],{x,xzero,xstar}
,Axes->False,Frame->True
,FrameLabel->{xlabel,label}
,LabelStyle->labelSize
,Exclusions->{xtilde}
,AspectRatio->1
,Epilog->If[OptionValue[ShowGuides],{
{Gray,Dashed,Line[{{xtilde,-(10^4)},{xtilde,+(10^4)}}]}
,{Text["ψ < 1",Scaled@{(xtilde/xstar),0.5},{0,-2},{0,1}]}
,{Text["ψ = 1",Scaled@{(xtilde/xstar),0.5},{0,2},{0,1}]}
},{}]
,Filling->0
,ImageSize->300
]
PlotInfo[{xzero_,xtilde_,xstar_},par:{a_,ab_,σ_,r_,ρ_}]:=Module[{α,β,γ},
α=αFun[par];β=βFun[par];γ=γFun[par];
Panel[
Grid[(({{"a", "=", ToString[a], "  ", "(OverBar[a])", "=", ToString[ab]}, {"σ", "=", ToString[σ], "  ", "", "", ""}, {"r", "=", ToString[r], "  ", "ρ", "=", ToString[ρ]}, {"(SuperStar[x])", "=", ToString[xstar], "  ", "(OverTilde[x])", "=", ToString[xtilde]}, {"α", "=", ToString[α], "  ", "β", "=", ToString[β]}, {"γ", "=", ToString[γ], "  ", "", "", ""}}))
,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}, {Center, Center, Left, Center, Center, Center, Left}}))
,Spacings->{0.2,0.8}
]]
]
PlotNumRule={0->"Off",1->qlabel,2->qplabel,3->Wlabel,4->Wplabel,5->flabel,6->Flabel,7->ψlabel,8->μxlabel,9->σxlabel,10->μqlabel,11->σqlabel};
Options[PlotBasicChoice]={ShowGuides->False};
PlotBasicChoice[plotnum_,{qn_,qpn_,Wn_,Wpn_,fn_,Fn_,ψn_,μxn_,σxn_,μqn_,σqn_},ran_,par:{a_,ab_,σ_,r_,ρ_},OptionsPattern[]]:=Switch[plotnum
,0,Graphics@Text[Style["Off",labelSize+smallBoost],{0,0}]
,1,PlotBasic[{qn,qlabel},ran,ShowGuides->OptionValue[ShowGuides]]
,2,PlotBasic[{qpn,qplabel},ran,ShowGuides->OptionValue[ShowGuides]]
,3,PlotBasic[{Wn,Wlabel},ran,ShowGuides->OptionValue[ShowGuides]]
,4,PlotBasic[{Wpn,Wplabel},ran,ShowGuides->OptionValue[ShowGuides]]
,5,PlotBasicFill[{fn,flabel},ran,ShowGuides->OptionValue[ShowGuides]]
,6,PlotBasicFill[{Fn,Flabel},ran,ShowGuides->OptionValue[ShowGuides]]
,7,PlotBasic[{ψn,ψlabel},ran,ShowGuides->OptionValue[ShowGuides]]
,8,PlotBasic[{μxn,μxlabel},ran,ShowGuides->OptionValue[ShowGuides]]
,9,PlotBasic[{σxn,σxlabel},ran,ShowGuides->OptionValue[ShowGuides]]
,10,PlotBasic[{μqn,μqlabel},ran,ShowGuides->OptionValue[ShowGuides]]
,11,PlotBasic[{σqn,σqlabel},ran,ShowGuides->OptionValue[ShowGuides]]
]

Solution and model output

Now that all the solver functions are defined, we can plug in some actual parameter values. As a baseline, we use the parameter values given in the text:

parn={0.11,0.07,0.1,0.05,0.06};
parnRule=MapThread[Rule,{par,parn}]

Automatic, interactive solver

pvOld=parn;
pvFirst=True;
Off[InterpolatingFunction::dmval]
Off[FindRoot::lstol]
Manipulate[
Catch[Module[{pv,ran,funs,plotax,plotbx,fail=False,epi},
pv={am,am abOveram,σm,rm,rm+ρMinusrm};
If[pvFirst||pv≠pvOld
,{okn,{qn,qpn,Wn,Wpn},{xzeron,xtilden,xstarn}}=Check[
SolveQWAuto[pv],
Throw["Unexpected message generated by SolveQWAuto","IKMFail"]
]
;{fn,Fn,{xf0n,xf1n}}=Check[
SolveErgodic[{qn,qpn,Wn,Wpn},{xtilden,xstarn},pv],
Throw["Unexpected message generated by SolveQWAuto","IKMFail"]
]
;{qppn,ψn,μxn,σxn,μqn,σqn}=SolveOthers[{qn,qpn,Wn,Wpn},{xtilden,xstarn},pv]
;pvOld=pv;pvFirst=False
];
funs={qn,qpn,Wn,Wpn,fn,Fn,ψn,μxn,σxn,μqn,σqn};
ran={xzeron,xtilden,xstarn};
epi=If[showi,Inset[PlotInfo[ran,pv],Scaled[{0.5,0}],{0,-1}],{}];
If[plota==0,plotax=plotb;plotbx=0,plotax=plota;plotbx=plotb];
Check[
If[plotbx≠0
,Show[GraphicsGrid[{{
PlotBasicChoice[plotax,funs,ran,pv,ShowGuides->showg],PlotBasicChoice[plotbx,funs,ran,pv,ShowGuides->showg]
}}],Epilog->epi]
,Show[PlotBasicChoice[plotax,funs,ran,pv,ShowGuides->showg],Epilog->epi]
],
Throw["Unexpected message generated by the plot function","IKMFail"]
]
]
,"IKMFail",{val,tag}|->Graphics@Text[Style["Solver failed: "<>val,labelSize+smallBoost],{0,0}]
]
,Style["Plot function selection",Bold]
,{{plota,1,"Plot A"},PlotNumRule}
,{{plotb,3,"Plot B"},PlotNumRule}
,{{showg,False,"Guides?"},{True,False}}
,{{showi,False,"Info?"},{True,False}}
,Delimiter
,Style["Model parameters",Bold]
,{{am,a/.parnRule,"a"},0.01,0.5}
,{{abOveram,ab/a/.parnRule,"(OverBar[a])/a"},0.01,0.99}
,{{σm,σ/.parnRule,"σ"},0.01,0.5}
,{{rm,r/.parnRule,"r"},0.01,0.15}
,{{ρMinusrm,ρ-r/.parnRule,"ρ - r"},0.001,0.15}
,ContinuousAction->False
,SaveDefinitions->True
]

Note: Root finder tolerance warning is turned off -- this may result in boundary conditions not being perfectly satisfied. This is appears to be a rare event, but if certainty about the validity of the solutions is needed, the boundary conditions should be checked by hand (or by plotting q' and W').