JMF simulations

Redistribution Recession sequels

Utilities

returns the eigenvector associated with the smallest (ie., most negative or least positive) eigenvalue.

eigenvectorsmall[m_]:=Last@First@Sort@Transpose@Eigensystem@m

rdiffpwr={(((-1+x_))^(y_))->(((-1))^y)(((1-x))^y),((((1/(-1+x_))))^(y_))->(((-1))^y)((((1/(1-x))))^y),(((-1))^(2 __))->1,(((-1))^(-2 __))->1};

distributepowers[expr_]:=FullSimplify[expr,β>0&&η>0&&γm>0&&γn>0&&n>0]//.(((a_ b_))^(c_))->((a^c)(b^c))//.rdiffpwr

font14[expr_]:=Style[expr,FontSize->14]

maxarginpiecewise[expr_/;Head@expr==Piecewise]:=Max@Cases[List@@#&/@Extract[#,Most/@Position[#,Alternatives@@{Equal,Less,LessEqual,Greater,GreaterEqual,Unequal}]]&@(Last/@First@expr)//Flatten,
_Integer|_Real]

NGM Analytical Setup

k[t]=date t capital
a=production function parameter (constant)
f[k[t],a]=date t output
c[t]=date t consumption
δ=date t depreciation rate (constant)
r[t]=date t capital rental rate
ρ=date t time preference rate (constant)
u[c[t],n[t]]=date t utility flow

Functional forms

f[k_,n_,a_]:=a (k^(1-α))(n^α)

u[c_,n_,m_]:=(σ/(σ-1))(c^((σ-1)/σ))-γn(η/(η+1))(((n-β m))^((η+1)/η))-γm(η/(η+1))(m^((η+1)/η))

Utility function when m is optimized

The composition of labor is constant.

First@Solve[Derivative[0,0,1][u][c,n,m]==0,m]

interpretγ=γ->γm γn ((((1/(((γm)^η)+(β^(1+η)) ((γn)^η)))))^(1/η));

u[c_,n_]:=(σ/(σ-1))(c^((σ-1)/σ))-γ(η/(η+1))(n^((η+1)/η))

u[c,n]==u[c,n,m]/.interpretγ/.First@Solve[Derivative[0,0,1][u][c,n,m]==0,m]//distributepowers

Equations

eqresourceconstraint=k'[t]==f[k[t],n[t],a]-c[t]-δ k[t];

eqintertemporalfoc=c'[t]==FullSimplify[Derivative[1,0][u][c[t],n[t]](( Derivative[1,0,0][f][k[t],n[t],a]-δ-ρ)/(-Derivative[2,0][u][c[t],n[t]]))]

eqlaborfoc=-((Derivative[0,1][u][c[t],n[t]])/(Derivative[1,0][u][c[t],n[t]]))==(1-τ[t]) Derivative[0,1,0][f][k[t],n[t],a]//FullSimplify

Inequalities

iqpositive=a>0&&γ>0&&η>0&&σ>0&&η>0&&c[t]>0&&n[t]>0&&k[t]>0&&0≤τ[t]<1 &&0<α<1;

Replacement rules

nofck[c_,k_,τ_]:=(((((1-τ)/γ)((α (k^(1-α)))/(c^(1/σ)))a))^(η/(1+(1-α)η)))

eqlaborfoc/.n[t]->nofck[c[t],k[t],τ[t]]//PowerExpand//distributepowers

rssnotation={τ[t]->τ,c[t]->css,k[t]->kss,n[t]->nss,r[t]->rss,τ'[t]->0,c'[t]->0,n'[t]->0,k'[t]->0,r'[t]->0};

Steady state

Steady-state version of equations and inequalities

iqsspositive=iqpositive/.rssnotation;

eqssintertemporalfoc=FullSimplify[eqintertemporalfoc/.rssnotation,iqsspositive]

eqsslaborfoc=FullSimplify[eqlaborfoc/.rssnotation,iqsspositive]

eqssresourceconstraint=eqresourceconstraint/.rssnotation

Steady-state behavior

rnoverk=FullSimplify@First@Quiet@Solve[eqssintertemporalfoc/.{css->(coverk kss),nss->(noverk kss)},noverk]//.rdiffpwr

rcoverk=FullSimplify@First@Solve[eqssresourceconstraint/.{css->(coverk kss),nss->(noverk kss)}/.rnoverk//PowerExpand,coverk]

rkss=First@Quiet@Solve[eqsslaborfoc/.{css->(coverk kss),nss->(noverk kss)}/.rnoverk/.rcoverk,kss]//PowerExpand//distributepowers//FullSimplify

Local behavior

localgradients=Grad[Last/@{eqintertemporalfoc,eqresourceconstraint}/.n[t]->nofck[c[t],k[t],τ],
{c[t],k[t]}]/.{c[t]->coverk kss,k[t]->kss}/.Flatten[{rcoverk,rkss}]//PowerExpand//distributepowers;

MatrixForm@localgradients

This depends only on α, δ, η, ρ, σ and not on a, γ, or τ

Initial impact on log labor

loglaborimpact=Dt[PowerExpand@Log@nofck[c,k,1-q],i]/.{Alternatives@@Dt[{a,α,γ,η,k,σ},i]->0}/.{Dt[q,i]->q Dt[lnq,i]}(* change in c beginning at c==css *)/.c->css/.Dt[css,i]->Dt[c,i]

Normalize γ and a so that kss == nss == 1 at a specific tax rate

rnorma=FullSimplify@First@Quiet@Solve[1==noverk kss/.rnoverk/.kss->1,a]

rnormγ[τ_]=distributepowers@PowerExpand@First@Quiet@Solve[1==kss /.rkss/.rnorma,γ]

Transition with constant parameters

Function that numerically solves for steady state and policy function

The tax rate in the calibration list must be the steady state tax rate

stationarysystem[calibration_List,leftbbd_,rightbbd_]:=If[Length@#<1,
Module[{rssnumerical=Inner[Rule,{css,kss,nss},({kss coverk,kss,kss noverk}//.Flatten[{rnoverk,rcoverk,rkss,calibration}]),List]},
{#(* return the calibration*),
Values@First@First@NDSolve[{csm'[k]==FullSimplify[Piecewise[{
{Apply[Divide,localgradients//.N[#]//eigenvectorsmall],
Abs[1-kss/k]<.001(* close to kss *)}},

(* otherwise: *)
((Last@eqintertemporalfoc)/(Last@eqresourceconstraint))]/.{k[t]->k,c[t]->csm[k],τ[t]->τ,n[t]->nofck[csm[k],k,τ]}//.#,
iqpositive/.{k[t]->k,c[t]->csm[k],τ[t]->τ}](* end diffeq *),
csm[kss]==css/.#(* bbd condn *)},
csm[k],{k,leftbbd kss/.#,rightbbd kss/.#}]
}&@Flatten@{rssnumerical,calibration}],
"Missing values for "<>ToString@#]&@Complement[{a,α,γ,δ,η,ρ,σ,τ},Keys@Flatten@calibration](* looking for missing params *)

stationarysystem[calibration_List]:=stationarysystem[calibration,0.7,1.5]

Plot example

Long run tax rate is 47%.  kss and nss are normalized to be the ss corresponding τ == 40%

Show@@Join[{Plot[Tooltip[Last@#,"Stable manifold"],{k,.9,1.1},PlotStyle->Black,AxesLabel->font14/@{"Capital","Consumption"}],
ListPlot[{Tooltip[{kss,css},"New steady state"]}/.First@#,PlotStyle->Black]}&@
stationarysystem@{#,τ->(47/100)},
List@ListPlot[{Tooltip[{kss,css},"Old steady state"]}/.First@stationarysystem@{#,τ->(40/100)},
PlotStyle->Blue]]&@{α->7/10,δ->1/20,η->1,σ->11/10,ρ->1/20,rnormγ[(40/100)],rnorma}

Transition with variable parameters

Utilities

kstart[kend_?NumericQ,taxpath_,calibration_List]:=Module[{tend=maxarginpiecewise@taxpath,
sm,rparams},
{rparams,sm}=stationarysystem@Join[{τ->(taxpath/.t->(2 tend))},calibration];
{k[t]/.#/.t->0,Join[#,{n[t]->nofck[c[t],k[t],τ[t]]/.#/.{τ[t]->taxpath}//.rparams,τ[t]->taxpath},Flatten@calibration]}&@First@NDSolve[{
eqintertemporalfoc,
eqresourceconstraint,
(* guess where the economy will be at time tend *)
k[tend]==kend,
(sm/.k->kend)==c[tend]}/.{n[t]->nofck[c[t],k[t],τ[t]]}/.{τ[t]->taxpath}//.rparams,
{c[t],k[t]},
{t,0,tend},
Method->{"DiscontinuityProcessing"->False}(* This option is necessary only if there are a number of parameter jumps *)]
]

timepaths[k0_,taxpath_,calibration_List]:=Module[{oneshot},
SetDelayed[oneshot[kend_],kstart[kend,taxpath,calibration]];
oneshot@NestWhile[(PrintTemporary@#;#-(((#1-k0)/(#2-#1))&[First@oneshot[#],First@oneshot[#+1/100]])(1/100))&,k0,Max[Abs[((#1-#3)/(#3))],Abs[((#2-#3)/(#3))]]>.001&,3,20]
]

processcalibration[descriptors_Association,calibration_List,endrates_List,numerics_:{False,False}]:=Module[{tmp=timepaths@@{1,Sequence@@calibration}},
Join[<|"timestamp"->Now|>,
descriptors,
<|"tastestechparams"->MapAt[#/.Take[Last@tmp,{5,-1}]&,Take[Last@tmp,{5,-1}],{All,2}],
"SS corresponding to time 0"->Take[First@stationarysystem@If[TrueQ@numerics[[1]],N,Identity]@Join[{τ->First@endrates},Take[Last@tmp,{5,-1}]],4],
"SS corresponding to LR"->Take[First@stationarysystem@If[TrueQ@numerics[[2]],N,Identity]@Join[{τ->Last@endrates},Take[Last@tmp,{5,-1}]],4],
"SM corresponding to LR"->Last@stationarysystem@If[TrueQ@numerics[[2]],N,Identity]@Join[{τ->Last@endrates},Take[Last@tmp,{5,-1}]],
"normalizing tax rate"->0.13,
"equilibrium time paths"->Inner[Rule,{c[t],k[t],n[t],τ[t]},{c[t],k[t],n[t],τ[t]}/.Last@tmp,List],
"tax rate plot"->calendarplot[τ[t]/.Last@tmp,{PlotRange->Full,PlotStyle->{Red,Thick}}]
|>]]

Plotting

timefunctiontocalendartime[timefunction_,startdate_List]:=timefunctiontocalendartime[timefunction,DateObject@startdate]

timefunctiontocalendartime[timefunction_,startdate_DateObject]:=
timefunctiontocalendartime[timefunction,startdate,120]

timefunctiontocalendartime[timefunction_,startdate_DateObject,gridsperyear_Integer]/;SameQ[Head@Head@timefunction,InterpolatingFunction]:=Module[{startgentime=Part[timefunction,0,1,1,1],
endgentime=Part[timefunction,0,1,1,-1],gentimes},
gentimes=Range[startgentime,endgentime,((endgentime-startgentime)/gridsperyear)];
TimeSeries[(timefunction/.t->#)&/@gentimes,{startdate+Quantity[gentimes,"Years"]}]]

timefunctiontocalendartime[timefunction_,startdate_DateObject,gridsperyear_Integer]:=Module[{startgentime=0,
endgentime=9,gentimes},
gentimes=Range[startgentime,endgentime,((endgentime-startgentime)/gridsperyear)];
TimeSeries[(timefunction/.t->#)&/@gentimes,{startdate+Quantity[gentimes,"Years"]}]]

calendarplot[timefunction_,options_:{}]:=DateListPlot[timefunctiontocalendartime[timefunction,{2007,12,1}],
Sequence@@options,ImageSize->500,Frame->{{True,False},{True,False}},
FrameTicks->{{#,DateString[#,{"MonthNameShort"," ","YearShort"}]}&/@(DateObject@{2007,12,1}+Quantity[1,"Months"]+Quantity[Range[0,9],"Years"]),Automatic}]

JMF final

multiplecalibrations5=Distribute[
{{Piecewise[{
{0.25+(t/2).12,t≥0&&t<2},
{0.3701,t≥2&&t<9}
},0.37],
Piecewise[{
{0.25+(t/2).12,t≥0&&t<2},
{0.37-((t-2)/(4 3)).12,t≥2&&t<6},
{0.37-(4/(4 3)).12+0.0001,t≥6&&t<9}
},0.37-(4/(4 3)).12],
Piecewise[{
{0.25+(t/2).12,t≥0&&t<2},
{0.37-((t-2)/4).12,t≥2&&t<6},
{0.2501,t≥6&&t<9}
},0.25]},

List@{α->7/10,δ->(65/1000)+(3/400)//Simplify,η->3/4,σ->135/100,ρ->1/200,nonprimesh->(1/10),rnormγ[(1/4)],rnorma}},List];

tmp5=processcalibration@@#&/@Transpose@{{<|"paperversion"->"JMF final",
"tastedescription"->"RR benchmark",
"techndescription"->"JMF benchmark",
"taxdescription"->"No recovery"|>,
<|"paperversion"->"JMF final",
"tastedescription"->"RR benchmark",
"techndescription"->"JMF benchmark",
"taxdescription"->"Partial recovery"|>,
<|"paperversion"->"JMF final",
"tastedescription"->"RR benchmark",
"techndescription"->"JMF benchmark",
"taxdescription"->"Full recovery"|>},
multiplecalibrations5,
{{(1/4),(37/100)},{(1/4),(1/4)*(1/3)+(2/3)*(37/100)},{(1/4),(1/4)}},
{{False,False},{False,False},{False,False}}};

Steady-state investment/adj-output ratio

multiplecalibrations5[[1,-1]]

{1-((css+(1-τ)α f[kss,nss,a](nonprimesh/(1-nonprimesh)))/(f[kss,nss,a]+(1-τ)α f[kss,nss,a](nonprimesh/(1-nonprimesh)))),1-(css/(f[kss,nss,a]))}//.{nss->nofck[coverk kss,kss,τ],css-> coverk kss}/.rcoverk/.rkss/.{τ->1/4}//.Flatten@multiplecalibrations5[[1,-1]]//N

write to file

LocalSymbol[NotebookDirectory[]<>"JMFsimulations"]=Join[LocalSymbol[$LocalSymbolDefaultDirectory<>"JMFsimulations"],
tmp5];

JMF Sensitivity analysis final

multiplecalibrations6={Piecewise[{
{0.25+(t/2).12,t≥0&&t<2},
{0.37-((t-2)/(4 3)).12,t≥2&&t<6},
{0.37-(4/(4 3)).12+0.0001,t≥6&&t<9}
},0.37-(4/(4 3)).12],#}&/@
{{α->7/10,δ->(65/1000)+(3/400)//Simplify,η->3/4,σ->135/100,ρ->1/200,nonprimesh->(1/10),rnormγ[(1/4)],rnorma},
{α->7/10,δ->(65/1000)+(3/400)//Simplify,η->3/8,σ->135/200,ρ->1/200,nonprimesh->(1/10),rnormγ[(1/4)],rnorma},
{α->7/10,δ->(65/1000)+(3/400)+(1/200)//Simplify,η->3/4,σ->135/100,ρ->1/200-((7/10)/200)//Simplify,nonprimesh->(1/10),rnormγ[(1/4)],rnorma},
{α->7/10,δ->(65/1000)+(3/400)-(1/200)//Simplify,η->3/4,σ->135/100,ρ->1/200+((7/10)/200)//Simplify,nonprimesh->(1/10),rnormγ[(1/4)],rnorma}};

Apply[SameQ,((α δ+ρ)/(1-α))/.Flatten/@multiplecalibrations6[[All,-1]]]

tmp6=processcalibration@@#&/@Transpose@{{
<|"paperversion"->"JMF final",
"tastedescription"->"RR benchmark",
"techndescription"->"JMF benchmark",
"taxdescription"->"Partial recovery"|>,
<|"paperversion"->"JMF final",
"tastedescription"->"smaller elasticities",
"techndescription"->"JMF benchmark",
"taxdescription"->"Partial recovery"|>,
<|"paperversion"->"JMF final",
"tastedescription"->"RR benchmark",
"techndescription"->"high TFP gr",
"taxdescription"->"Partial recovery"|>,
<|"paperversion"->"JMF final",
"tastedescription"->"RR benchmark",
"techndescription"->"zero TFP gr",
"taxdescription"->"Partial recovery"|>},
multiplecalibrations6,
{{(1/4),(1/4)*(1/3)+(2/3)*(37/100)},{(1/4),(1/4)*(1/3)+(2/3)*(37/100)},{(1/4),(1/4)*(1/3)+(2/3)*(37/100)},{(1/4),(1/4)*(1/3)+(2/3)*(37/100)}},
{{False,False},{False,False},{False,False},{False,False}}};

write to file

LocalSymbol[NotebookDirectory[]<>"JMFsimulations"]=Join[LocalSymbol[$LocalSymbolDefaultDirectory<>"JMFsimulations"],
Rest@tmp6];