(*This file will compute the function m_0 as given in Figures 1 and 2 from the paper "Engodenous Current Coupons"*)

(*CIR Parameters*)
κ=0.25;
θ=.06;
σ=.1;

(*invariant density for the CIR*)
m[x_]:= ((((2×(κ/(σ^2))))^(2×κ×(θ/(σ^2))))/(Gamma[2×κ×(θ/(σ^2))]))×(x^(2×κ×(θ/(σ^2))-1))×Exp[-2×(κ/(σ^2))×x];

(*Bond Pricing formula:
	f(t,x) = E^x[ e^{-\int_0^t r_u du}] = e^{-xCC(t) - AA(t)} *)
tgamma = (1/2)×((((κ^2)+2×(σ^2)))^(1/2));
CC[t_]:=((Sinh[tgamma×t])/(tgamma×Cosh[tgamma×t]+(1/2)×κ×Sinh[tgamma×t]));
AA[t_]:=((-2×κ×θ)/(σ^2))×Log[tgamma×((Exp[(1/2)×κ×t])/(tgamma×Cosh[tgamma×t]+(1/2)×κ×Sinh[tgamma×t]))];
f[t_,x_]:=Exp[-x×CC[t]-AA[t]];

(* function for finding m0 when γ_0(x)= γ > 0.  We wish to solve g[m,x,T,γ] = 0 to obtain a function m = m(x) = m(x; T,γ)*)g[m_,x_,T_,γ_]:=NIntegrate[Exp[-γ×t]×f[t,x](1+γ×((1-Exp[-m×(T-t)])/m)),{t,0,T}]-((1-Exp[-m×T])/m);
m0[x_,γ_]:=m/.FindRoot[g[m,x,30,γ],{m,x}][[1]]; 
(*hardcode 30 years since this works better numerically*)

(* function for finding m0 when γ_0(x) = 0. Here, we have an explicit formula using the ProductLog function.*)
hatg[x_,T_]:=(1/T)×NIntegrate[f[t,x],{t,0,T}];
hatm0[x_,T_]:=(1/T)×((1+hatg[x,T]×ProductLog[((-1)/(hatg[x,T]))×Exp[((-1)/(hatg[x,T]))]])/(hatg[x,T]));

(*Turn off NIntegrate and FindRoot warnings...*)
Off[NIntegrate::inumr];
Off[NIntegrate::slwcon];
Off[NIntegrate::ncvb];
Off[FindRoot::lstol];

(*Obtain m_0(x) when γ_0(x) = 0 and γ_0(x)=0.045*)
(*table outputs the pair (x,m_0(x)) along a mesh with width 1/1000 for 0≤x≤.5*)
gam0Table = Table[{i/1000,hatm0[i/1000,30]},{i,1,500}];
gam045Table = Table[{i/1000,m0[i/1000,0.045]},{i,1,500}];

(*the following two tables strip out just m_0(x) from the above tables.*)gam0AltTable = Table[gam0Table[[i,2]],{i,1,500}];
gam045AltTable = Table[gam045Table[[i,2]],{i,1,500}];

(*Plot the functions*)
ListLinePlot[{gam0Table,gam08Table},PlotStyle->{Dashing[None],Dashing[Tiny]}]