Pension Plans with Mathematica And Python

Seth J. Chandler 2013

Wolfram Technology Conference 2013

Using Mathematica to generate benefit side

Constants and functions

Needs["Pensions`"]

baseDirectory=(*"/Users/sethchandler/Dropbox/Novim Pensions";*)"/Users/sethchandler/Documents/TemporaryPensions/";

commonDirectory=(*"/Users/sethchandler/Dropbox/Pensions/3 - Modeling";*)"/Users/sethchandler/Documents/TemporaryPensions/";

notWorkingAnymoreHazard[survival_,survivalVectorInitialAge_,separationHazardRaggedRight_,initialEntryAge_,entryAge_]:=balancedMultipleDecrementHazardVectorC[survivalToHazardC[truncateSurvivalVectorC3[survival,survivalVectorInitialAge,entryAge]],separationHazardRaggedRight[[1+entryAge-initialEntryAge]]]

AntiDiagonalUpperTriangularToUpperTriangular=Compile[{{m,_Real,2}},MapIndexed[{row,index}|->With[{i=index[[1]]-1},RotateRight[row,i]],m],CompilationTarget->"C",CompilationOptions->{"InlineCompiledFunctions"->True,"ExpressionOptimization"->True},RuntimeOptions->"Speed"]

UpperTriangularToAntiDiagonalUpperTriangular=Compile[{{m,_Real,2}},MapIndexed[{row,index}|->With[{i=index[[1]]-1},RotateLeft[row,i]],m],CompilationTarget->"C",CompilationOptions->{"InlineCompiledFunctions"->True,"ExpressionOptimization"->True},RuntimeOptions->"Speed"]

workforceEvolveStep=Compile[{{persistence,_Real,2},{workforceNow,_Real,2},{distribution,_Real,1},{growth,_Real}},
Module[{totalNow=Total[Flatten[workforceNow]],newBeforeRotation,totalNext,newHires},
(*Assert[Total[distribution]]==1;*)
newBeforeRotation=persistence*workforceNow;
totalNext=Total[Flatten[newBeforeRotation]];
newHires=((1+growth)totalNow-totalNext)*distribution;
DiagonalMatrix[newHires]+Map[row|->RotateRight[row],persistence*workforceNow]
],CompilationTarget->"C",CompilationOptions->{"InlineCompiledFunctions"->True,"ExpressionOptimization"->True},RuntimeOptions->"Speed"
];

workforceEvolve=Compile[{{persistence,_Real,2},{workforceNow,_Real,2},{distribution,_Real,1},{growth,_Real},{steps,_Integer}},
NestList[workforce|->workforceEvolveStep[persistence,workforce,distribution,growth],workforceNow,steps],CompilationTarget->"C",
CompilationOptions->{"InlineCompiledFunctions"->True,"ExpressionOptimization"->True},RuntimeOptions->"Speed"
];

Data entry

entryAge=25;

serviceDomain={0,55};

retirementDomain={20,75};

survivalTableInitialAge=20;

killOffAge=120;

discountRate=0.05;

normalCostDiscountRate=0.05 (* check with Josh on this *);

COLARate=0.;

inflation=0.025;

fasBasisYears=5;

employeeContributionRatesList={{0,1000000,0.064}};

employerContributionRatesList={{0,1000000,0}};

ϕ=1.;

survival={1.`,0.999878`,0.9997590145180001`,0.9996390434362579`,0.9995140885558282`,0.9993821526961388`,0.9992432385769141`,0.9990903543614119`,0.9989314989950684`,0.9987636785032372`,0.9985868973321421`,0.9983881785395731`,0.9981435734358308`,0.9978640932352688`,0.9975587468227388`,0.997231547553781`,0.9968845109752322`,0.9965186543597043`,0.9961359911964301`,0.995734548391978`,0.9953123569434598`,0.9948525226345519`,0.9943481324055762`,0.9937932861476939`,0.9931830970699992`,0.9925136916625741`,0.9917980892908854`,0.9909897738481134`,0.9900423876243146`,0.9889087890904849`,0.987543106052751`,0.9858998343242792`,0.9841439467193476`,0.982178611257749`,0.9799451370957488`,0.977388460233066`,0.9744494531331452`,0.9710583690362419`,0.967206180486275`,0.9629001785707502`,0.9580827889773608`,0.9527089026139868`,0.9467458975925258`,0.9401603331288723`,0.9329239190447793`,0.9249949986568178`,0.91632317054441`,0.9068630501317095`,0.8965755956910154`,0.8854069534954924`,0.8732697949769761`,0.8600441239320499`,0.8458774771226412`,0.8303793099868002`,0.8138057693387737`,0.7958035719152308`,0.7767560134214397`,0.7562721805915029`,0.7338600545196737`,0.7098870481186795`,0.6842842618412313`,0.6570045854586688`,0.6280149151298905`,0.5972949375413968`,0.5648540575987109`,0.5307340482494607`,0.4942917255604597`,0.45566925300034206`,0.41515296570056365`,0.3740374614364769`,0.33206858807599704`,0.2908914190173973`,0.25128044359837926`,0.21398389247684624`,0.17894531398711752`,0.14751822216344398`,0.11995075439665039`,0.09628231144061213`,0.07596905450211754`,0.05929308734835772`,0.045845237258488145`,0.035174162298490665`,0.026562331446302`,0.019802271217881032`,0.014533995773991093`,0.010478211583280007`,0.0074068801168353056`,0.005126960941192113`,0.0034723624714458885`,0.0023006450067187163`,0.0014918670584767918`,0.0009479069672161595`,0.0005912607624289482`,0.0003630252392199378`,0.00022017226641021773`,0.00013247589132089674`,0.00007948553479253804`,0.00004769132087552282`,0.000028614792525313693`,0.000017168875515188216`,0.00001030132530911293`};(*Import["/Users/sethchandler/Dropbox/pensions/3 - Modeling/TRS Pension Wealth Accrual (JM 04162013).xlsx",{"Sheets","Female Actuarial Table"}][[3;;-1,3]]*)

multipleHazardTX=Module[{TXseparationTableRaw=Import[FileNameJoin[{commonDirectory,"TRS Dec Table.xlsx"}],{"Sheets","TRS Decrement Table"}],TXseparationTableRaw2,TXseparationTableRaw3,TXseparationTableRaw4},
TXseparationTableRaw2=Drop[TXseparationTableRaw,1,1];
TXseparationTableRaw3=Join[Table[TXseparationTableRaw2[[1]],{5}],TXseparationTableRaw2];
TXseparationTableRaw4=MapIndexed[PadRight[#,56+1-#2,1.]&,Replace[TXseparationTableRaw3,{a___?NumericQ,b_?(Not[NumericQ[#]]&),c___?(Not[NumericQ[#]]&)}:>{a,1},{1}]];
Table[ReplacePart[PadLeft[notWorkingAnymoreHazard[survival,20,TXseparationTableRaw4,20,entryAge],56],-1->1],{entryAge,20,75}]
];

MatrixPlot[multipleHazardTX,PlotLegends->Automatic]

txReplRateImage=;

txReplRate=ImageData[txReplRateImage];

salaryScale={44130.`,44130.`,44130.`,44233.`,44336.`,44439.`,44542.`,44748.`,44748.`,45057.`,45366.`,45675.`,46293.`,46911.`,47529.`,48147.`,48765.`,49383.`,50001.`,50619.`,51237.`,51855.`,52473.`,53091.`,53709.`,54327.`,54945.`,55563.`,56181.`,57417.`,58035.`,58653.`,59889.`,59889.`,60507.`,61125.`,62361.`,62919.`,64833.`,64833.`,64833.`,64833.`,64833.`,64833.`,64833.`,64833.`,64833.`,64833.`,64833.`,64833.`,64833.`,64833.`,64833.`,64833.`,64833.`,64833.`};(*Import["/Users/sethchandler/Dropbox/pensions/3 - Modeling/TRS Pension Wealth Accrual (JM 04162013).xlsx",{"Sheets","sal grid"}][[2;;-1,2]]*)

wfImage=;

workforceSnapshot=ImageData[wfImage];

workforceReplacementDistribution=With[{hiresThisYear=Diagonal[workforceSnapshot]},N[(hiresThisYear/(Total[hiresThisYear]))]]

Additional Functionality

bc[salaryScale_,fasBasisYears_,survival_,survivalTableInitialAge_,discountRate_,inflation_,replacementRateTable_,basedExclusivelyOnRetirementQ_,compoundQ_,COLARate_,lastSomethingWithoutIncrement_,COLACap_,killOffAge_,employeeContributionRatesList_,employerContributionRatesList_,serviceDomain_,retirementDomain_,person:{entryAge_,yos_}]:=Module[{age=entryAge+yos,mmcd=matrixMatrixCommonDomainC[{{serviceDomain,retirementDomain},{{0,retirementDomain[[2]]-entryAge},retirementDomain}}],eec,erc,αRR,αSDAF,αΩ,retirementAgeIndexStar,bim,z,salary,finalAverageSalary,cumulativeInvestedContributions,refundRight},
salary=Take[salaryVectorFromScaleEntryAgeAndYOSC[salaryScale,inflation,{entryAge,yos}],1+retirementDomain[[2]]-entryAge];
eec=contributionsC[salary,employeeContributionRatesList];
erc=contributionsC[salary,employerContributionRatesList];
αRR=trimMatrixByCommonDomainC[replacementRateTable,{serviceDomain,retirementDomain},mmcd];
αSDAF=trimMatrixByCommonDomainC[sDAFMatrixC[afRetirementAgeVectorC[survival,survivalTableInitialAge,discountRate,colaMatrixC[basedExclusivelyOnRetirementQ,compoundQ,COLARate,lastSomethingWithoutIncrement,COLACap,retirementDomain,killOffAge],retirementDomain],survivalDiscountRetirementAgeMatrixC[survival,survivalTableInitialAge,discountRate,{entryAge,retirementDomain[[2]]},retirementDomain]],{{0,retirementDomain[[2]]-entryAge},retirementDomain},mmcd];

αΩ=omegaC[αRR,αSDAF];
retirementAgeIndexStar=bestPositionsInMatrixByRowC@αΩ;
z=retirementRRSDAFΩStar[αRR,αSDAF,retirementDomain];

finalAverageSalary=fasVectorC[salary,fasBasisYears];
cumulativeInvestedContributions=Most@cumulativeInvestedFlowsC[eec,discountRate];
refundRight=Most@refundRightC[eec,erc,1,discountRate];
benefitSideC[z,salary,finalAverageSalary,cumulativeInvestedContributions,refundRight,inflation]
]

bothSides[salaryScale_,fasBasisYears_,survival_,survivalTableInitialAge_,discountRate_,inflation_,replacementRateTable_,basedExclusivelyOnRetirementQ_,compoundQ_,COLARate_,lastSomethingWithoutIncrement_,COLACap_,killOffAge_,employeeContributionRatesList_,employerContributionRatesList_,serviceDomain_,retirementDomain_,person:{entryAge_,yos_},notWorkingHazard_,discountRateForEmployer_,paymentTimeRelativeToHazardSalary_:1,paymentTimeRelativeToHazardPension_:1]:=Module[{age=entryAge+yos,mmcd=matrixMatrixCommonDomainC[{{serviceDomain,retirementDomain},{{0,retirementDomain[[2]]-entryAge},retirementDomain}}],eec,erc,αRR,αSDAF,αΩ,retirementAgeIndexStar,bim,z,salary,finalAverageSalary,cumulativeInvestedContributions,refundRight},
salary=Take[salaryVectorFromScaleEntryAgeAndYOSC[salaryScale,inflation,{entryAge,yos}],1+retirementDomain[[2]]-entryAge];
eec=contributionsC[salary,employeeContributionRatesList];
erc=contributionsC[salary,employerContributionRatesList];
αRR=trimMatrixByCommonDomainC[replacementRateTable,{serviceDomain,retirementDomain},mmcd];
αSDAF=trimMatrixByCommonDomainC[sDAFMatrixC[afRetirementAgeVectorC[survival,survivalTableInitialAge,discountRate,colaMatrixC[basedExclusivelyOnRetirementQ,compoundQ,COLARate,lastSomethingWithoutIncrement,COLACap,retirementDomain,killOffAge],retirementDomain],survivalDiscountRetirementAgeMatrixC[survival,survivalTableInitialAge,discountRate,{entryAge,retirementDomain[[2]]},retirementDomain]],{{0,retirementDomain[[2]]-entryAge},retirementDomain},mmcd];

αΩ=omegaC[αRR,αSDAF];
retirementAgeIndexStar=bestPositionsInMatrixByRowC@αΩ;
z=retirementRRSDAFΩStar[αRR,αSDAF,retirementDomain];

finalAverageSalary=fasVectorC[salary,fasBasisYears];
cumulativeInvestedContributions=Most@cumulativeInvestedFlowsC[eec,discountRate];
refundRight=Most@refundRightC[eec,erc,1,discountRate];
benefitAndCostSideSingleStepC[z,salary,finalAverageSalary,cumulativeInvestedContributions,refundRight,inflation,notWorkingHazard,discountRateForEmployer,paymentTimeRelativeToHazardSalary,paymentTimeRelativeToHazardPension]
]

simpleUserInterfaceBothSides[colaRate_,deflateForInflation_,discountRate_,employeeContributionRate_,entryAge_,fasYears_,inflation_,showTooltips_]:=Module[{data$=bothSides[salaryScale,fasYears,survival,survivalTableInitialAge,discountRate,inflation,txReplRate,True,True,colaRate,0.`,1000,killOffAge,{{0,1000000,employeeContributionRate}},employerContributionRatesList,serviceDomain,retirementDomain,{entryAge,0},multipleHazardTX[[entryAge-19]],discountRate,1,1],tips$},tips$=Transpose[data$[[{1,2,3,4,5,8,9,10,-2,-1}]]];Labeled[ListLinePlot[If[showTooltips,MapThread[Function[{datum,tip},Tooltip[datum,Grid[Transpose[{{"salary","final average salary","cumulative invested contributions","refund right","optimal retirement age","optimal replacement rate","first year payment","present  value of annuity","normal cost percentage","normal cost"},tip}],Dividers->All,BaseStyle->{FontSize->11,LineIndent->0},ItemSize->15]]],{If[deflateForInflation,data$[[{-5,-4}]],data$[[{-7,-6}]]],{tips$,tips$}},2],If[deflateForInflation,data$[[{-5,-4}]],data$[[{-7,-6}]]]],Axes->False,Frame->True,FrameLabel->{"years of service","$"},PlotLegends->{Style[If[deflateForInflation,"deflated ",""]<>"accrued pension wealth",{FontSize->14,FontFamily->"Swiss"}],Style[If[deflateForInflation,"deflated ",""]<>"accrued net pension wealth",{FontSize->14,FontFamily->"Swiss"}]},PlotStyle->Thick,BaseStyle->{FontSize->14,FontFamily->"Swiss"},PlotMarkers->{Automatic,Tiny},ImageSize->500,ImagePadding->{{80,10},{50,10}}],Row[{"Normal cost percentage is ",data$[[-2,1]]},BaseStyle->{Background->Lighter@ColorData[1][3],FontSize->14,FontFamily->"Swiss"}]]]

Manipulate[simpleUserInterfaceBothSides[colaRate,deflateForInflation,discountRate,employeeContributionRate,entryAge,fasYears,inflation,showTooltips],
{{showTooltips,False,"provide data tooltips"},{True,False}},
{{deflateForInflation,True,"show deflated values to take account of inflation"},{True,False}},
Row[{"initial pension annuity is replacement rate multiplied by the average salary for the final ... ",Control@{{fasYears,5,""},Range[7],ControlType->SetterBar},"... years"}],
{{discountRate,0.05,"discount rate"},0.01,0.1,0.005,Appearance->"Labeled"},
{{inflation,0.025,"inflation rate"},0.,0.1,0.005,Appearance->"Labeled"},
{{entryAge,25,"entry age"},20,60,1,Appearance->"Labeled"},
{{colaRate,0,"COLA Rate"},0,0.1,0.005,Appearance->"Labeled"},
{{employeeContributionRate,0.064,"employee contribution rate"},0.0,0.08,0.001,Appearance->"Labeled"}]

DynamicModule[{colaRate=0,deflateForInflation=True,discountRate=0.055`,employeeContributionRate=0.064`,entryAge=29,fasYears=3,inflation=0.025`,showTooltips=True},simpleUserInterfaceBothSides[colaRate,deflateForInflation,discountRate,employeeContributionRate,entryAge,fasYears,inflation,showTooltips]]

Graphics produced for video

Labeled[ListLinePlot[With[{t=Table[bothSides[salaryScale,fasBasisYears,survival,survivalTableInitialAge,discountRate,inflation,txReplRate,True,True,0.,0.,1000,killOffAge,employeeContributionRatesList,employerContributionRatesList,serviceDomain,retirementDomain,{ea,0},multipleHazardTX[[ea-19]],0.05,1,1],{ea,20,75,5}]},
Map[Thread[List[Range[1+75-Length[#],75],#]]&,Part[t,All,-4]]],Axes->False,Frame->True,FrameLabel->{"age","value"},ImageSize->600,PlotRange->All,Joined->True,PlotLegends->Range[20,75,5],ImageSize->400,PlotStyle->Thick,BaseStyle->{FontFamily->"Swiss",FontSize->16}],
Style["Deflated Net Pension Wealth Trajectories\nAs A Function of Entry Age",{FontFamily->"Swiss",FontSize->16}]
]

Labeled[ListLinePlot[With[{t=Table[bothSides[salaryScale,fasBasisYears,survival,survivalTableInitialAge,discountRate,inflation,txReplRate,True,True,0.,0.,1000,killOffAge,employeeContributionRatesList,employerContributionRatesList,serviceDomain,retirementDomain,{ea,0},multipleHazardTX[[ea-19]],0.05,1,1],{ea,20,75,5}]},
Map[Thread[List[Range[1+75-Length[#],75],#]]&,Part[t,All,-1]]],Axes->False,Frame->True,FrameLabel->{"age","value"},ImageSize->600,PlotRange->All,Joined->True,PlotLegends->Range[20,75,5],ImageSize->400,PlotStyle->Thick,BaseStyle->{FontFamily->"Swiss",FontSize->16}],
Style["Normal Cost Trajectories\nAs A Function of Entry Age",{FontFamily->"Swiss",FontSize->16}]
]

Manipulate[With[{t=Table[bothSides[salaryScale,fasBasisYears,survival,survivalTableInitialAge,discountRate,inflation,txReplRate,True,True,0.,0.,1000,killOffAge,employeeContributionRatesList,employerContributionRatesList,serviceDomain,retirementDomain,{ea,0},multipleHazardTX[[ea-19]],0.05,1,1],{ea,20,75,1}],labels={"salary","FAS","cumulative invested contributions","refund right","retirement age","Ωstar","SDAFStar","Replacement Rate Star","first year annuity payments star","pv annuity as of separation star","annuityQ","pension wealth","net pension wealth","pension wealth deflated","net pension wealth deflated","discounted net pension wealth","normal cost %","normal cost"}},Labeled[f[MapIndexed[Tooltip[Thread[List[Range[1+75-Length[#],75],#]],Range[20,75,1][[#2[[1]]]]]&,Part[t,All,part]],Axes->False,Frame->True,FrameLabel->{"age","value"},ImageSize->600,PlotRange->All,Joined->True,PlotLegends->Range[20,75,1],ImageSize->400,BaseStyle->{FontFamily->"Swiss",FontSize->16}],labels[[part]]]
],{part,1,18,1},{f,{ListLinePlot,ListLogPlot}}]

Aggregate cost trajectories

Manipulate[Labeled[MatrixPlot[Table[PadRight[bothSides[salaryScale,fasYears,survival,survivalTableInitialAge,discountRate,inflation,txReplRate,True,True,colaRate,0.`,1000,killOffAge,{{0,1000000,employeeContributionRate}},employerContributionRatesList,serviceDomain,retirementDomain,{entryAge,0},multipleHazardTX[[entryAge-19]],discountRate,1,1][[-1]],56],{entryAge,20,75,1}],FrameTicks->{Table[{i,i+20},{i,0,55,5}],Table[{i+1,i},{i,0,55,5}]},FrameLabel->{"entry age","complete years of service"},Mesh->{Range[0,55,5],Range[0,55,5]},BaseStyle->{FontSize->12,FontFamily->"Swiss"}],Style["Normal cost by entry age and completed years of service",{FontSize->12,FontFamily->"Swiss"}]],Row[{"initial pension annuity is replacement rate multiplied by the average salary for the final ... ",Control@{{fasYears,5,""},Range[7],ControlType->SetterBar},"... years"}],
{{discountRate,0.05,"discount rate"},0.01,0.1,0.005,Appearance->"Labeled"},
{{inflation,0.025,"inflation rate"},0.,0.1,0.005,Appearance->"Labeled"},
{{colaRate,0,"COLA Rate"},0,0.1,0.005,Appearance->"Labeled"},
{{employeeContributionRate,0.064,"employee contribution rate"},0.0,0.08,0.001,Appearance->"Labeled"},SynchronousUpdating->False
]

Module[{colaRate=0,discountRate=0.05`,employeeContributionRate=0.064`,fasYears=5,inflation=0.025`},Labeled[MatrixPlot[AntiDiagonalUpperTriangularToUpperTriangular[Table[PadRight[bothSides[salaryScale,fasYears,survival,survivalTableInitialAge,discountRate,inflation,txReplRate,True,True,colaRate,0.`,1000,killOffAge,{{0,1000000,employeeContributionRate}},employerContributionRatesList,serviceDomain,retirementDomain,{entryAge,0},multipleHazardTX[[entryAge-19]],discountRate,1,1][[-1]],56],{entryAge,20,75,1}]],FrameTicks->{Table[{i,i+20},{i,0,55,5}],Table[{i,i+20},{i,0,55,5}]},FrameLabel->{"entry age","age"},Mesh->{Range[0,55,5],Range[0,55,5]},BaseStyle->{FontSize->12,FontFamily->"Swiss"}],Style["Normal cost by entry age and completed years of service",{FontSize->12,FontFamily->"Swiss"}]]
]

Module[{colaRate=0,discountRate=0.05`,employeeContributionRate=0.064`,fasYears=5,inflation=0.025`},Labeled[MatrixPlot[AntiDiagonalUpperTriangularToUpperTriangular[Table[PadRight[bothSides[salaryScale,fasYears,survival,survivalTableInitialAge,discountRate,inflation,txReplRate,True,True,colaRate,0.`,1000,killOffAge,{{0,1000000,employeeContributionRate}},employerContributionRatesList,serviceDomain,retirementDomain,{entryAge,0},multipleHazardTX[[entryAge-19]],discountRate,1,1][[1]],56],{entryAge,20,75,1}]],FrameTicks->{Table[{i,i+20},{i,0,55,5}],Table[{i,i+20},{i,0,55,5}]},FrameLabel->{"entry age","age"},Mesh->{Range[0,55,5],Range[0,55,5]},BaseStyle->{FontSize->12,FontFamily->"Swiss"}],Style["Salary by entry age and completed years of service",{FontSize->12,FontFamily->"Swiss"}]]
]

Timing[Grid[𝒞=Table[PadRight[bothSides[salaryScale,fasBasisYears,survival,survivalTableInitialAge,discountRate,inflation,txReplRate,True,True,0.,0.,1000,killOffAge,employeeContributionRatesList,employerContributionRatesList,serviceDomain,retirementDomain,{ea,0},multipleHazardTX[[ea-19]],discountRate,1,1][[-1]],56],{ea,20,75,1}],BaseStyle->6]];

Labeled[MatrixPlot[workforceSnapshot,FrameTicks->{Table[{i+1,i+20},{i,0,55,5}],Table[{i+1,i+20},{i,0,55,5}]},FrameLabel->{"entry age","age"},BaseStyle->{FontSize->12,FontFamily->"Swiss"}],Style["entry age and age distribution of Texas teachers",{FontSize->12,FontFamily->"Swiss"}]
]

Magnify[Row[{Labeled[MatrixPlot[workforceSnapshot,FrameTicks->None],Style["workforce",10]],"×",Labeled[MatrixPlot[UpperTriangularize@(1-multipleHazardTX),FrameTicks->None],Style["persistence",10]],"=",Labeled[MatrixPlot[workforceSnapshot×UpperTriangularize@(1-multipleHazardTX),FrameTicks->None],Style["new workforce\nbefore aging and replacement",10]]},Spacer[12]],1.5]

Magnify[Row[{Labeled[MatrixPlot[workforceSnapshot,FrameTicks->None],Style["workforce",10]],"×",Labeled[MatrixPlot[UpperTriangularize@(1-multipleHazardTX),FrameTicks->None],Style["persistence",10]],"=",Labeled[MatrixPlot[Transpose@ReplacePart[Transpose@Map[RotateRight,workforceSnapshot×UpperTriangularize@(1-multipleHazardTX)],1->ConstantArray[0,56]],FrameTicks->None,PlotLabel->"rotated right"],Style["new workforce\nafter aging before replacement",10]]},Spacer[12]],1.5]

ListLinePlot[(Total[Flatten[workforceSnapshot]]*1.03-Total[Flatten[Transpose@ReplacePart[Transpose@Map[RotateRight,workforceSnapshot×UpperTriangularize@(1-multipleHazardTX)],1->ConstantArray[0,56]]]])*((Diagonal[workforceSnapshot])/(Total[Diagonal[workforceSnapshot]])),DataRange->{20,75},PlotStyle->Thick,PlotMarkers->Automatic,Axes->False,Frame->True,FrameLabel->{"age","new workers"},BaseStyle->{FontSize->14,FontFamily->"Swiss"},ImageSize->700]

Magnify[Row[{Labeled[MatrixPlot[DiagonalMatrix[(Total[Flatten[workforceSnapshot]]*1.03-Total[Flatten[Transpose@ReplacePart[Transpose@Map[RotateRight,workforceSnapshot×UpperTriangularize@(1-multipleHazardTX)],1->ConstantArray[0,56]]]])*((Diagonal[workforceSnapshot])/(Total[Diagonal[workforceSnapshot]]))],FrameTicks->None],"new workers"],"+",Labeled[MatrixPlot[Transpose@ReplacePart[Transpose@Map[RotateRight,workforceSnapshot×UpperTriangularize@(1-multipleHazardTX)],1->ConstantArray[0,56]],FrameTicks->None,PlotLabel->"rotated right"],Style["new workforce\nafter aging before replacement",10]],"=",Labeled[MatrixPlot[DiagonalMatrix[(Total[Flatten[workforceSnapshot]]*1.03-Total[Flatten[Transpose@ReplacePart[Transpose@Map[RotateRight,workforceSnapshot×UpperTriangularize@(1-multipleHazardTX)],1->ConstantArray[0,56]]]])*((Diagonal[workforceSnapshot])/(Total[Diagonal[workforceSnapshot]]))]+Transpose@ReplacePart[Transpose@Map[RotateRight,workforceSnapshot×UpperTriangularize@(1-multipleHazardTX)],1->ConstantArray[0,56]],FrameTicks->None],"new workforce"]},Spacer[6]],1.5]

Labeled[MatrixPlot[UpperTriangularize@(1-multipleHazardTX),FrameTicks->{Table[{i+1,i+20},{i,0,55,5}],Table[{i+1,i+20},{i,0,55,5}]},FrameLabel->{"entry age","age"},BaseStyle->{FontSize->12,FontFamily->"Swiss"}],Style["persistence of Texas teachers",{FontSize->12,FontFamily->"Swiss"}]
]

Module[{tensor=workforceEvolve[1-multipleHazardTX,workforceSnapshot,((Diagonal[workforceSnapshot])/(Total[Diagonal[workforceSnapshot]])),0.03,20],max},max=Max[Flatten[tensor]];
Image3D[workforceTensor=tensor/max,Axes->True,AxesLabel->{"entry age","age","future"},ImageSize->600]]

𝒞=Table[PadRight[bothSides[salaryScale,fasBasisYears,survival,survivalTableInitialAge,discountRate,inflation,txReplRate,True,True,0.,0.,1000,killOffAge,employeeContributionRatesList,employerContributionRatesList,serviceDomain,retirementDomain,{ea,0},multipleHazardTX[[ea-19]],discountRate,1,1][[-1]],56],{ea,20,75,1}];

Functionality and Data related to aggregation

wh=workforceEvolve[1-multipleHazardTX,workforceSnapshot,((Diagonal[workforceSnapshot])/(Total[Diagonal[workforceSnapshot]])),0.0,50];

notAdjustedForInflationCostTensor=Compile[{{workforceHistory,_Real,3},{costMatrix,_Real,2}},MapIndexed[{workforce,index}|->With[{t=index[[1]]-1},UpperTriangularToAntiDiagonalUpperTriangular[workforce]*costMatrix],workforceHistory],CompilationTarget->"C",
CompilationOptions->{"InlineCompiledFunctions"->True,"ExpressionOptimization"->True},RuntimeOptions->"Speed"
];

generateTensorIndices=Compile[{{costMatrixDimensions,_Integer,1},{yearsToProject,_Integer}},Outer[List,Range@costMatrixDimensions[[1]],Range@costMatrixDimensions[[2]],Range@yearsToProject],CompilationTarget->"C",
CompilationOptions->{"InlineCompiledFunctions"->True,"ExpressionOptimization"->True},RuntimeOptions->"Speed"
];

adjustedForInflationCostTrajectory=Compile[{{unadjustedCostTensor,_Real,3},{inflation,_Real}},
Module[{dims=Dimensions[unadjustedCostTensor],indices},
indices=generateTensorIndices[Take[dims,-2],dims[[1]]];
Map[Total[Flatten[#]]&,Take[MapThread[#1*((1+inflation)^(#2[[3]]+#2[[1]]-#2[[2]]))&,{unadjustedCostTensor,indices},3],Length[unadjustedCostTensor]]
]
],
CompilationTarget->"C",
CompilationOptions->{"InlineCompiledFunctions"->True,"ExpressionOptimization"->True},RuntimeOptions->"Speed"
];

costTrajectory=Compile[{{workforceHistory,_Real,3},{costMatrix,_Real,2},{inflation,_Real}},
With[{notAdjusted=notAdjustedForInflationCostTensor[workforceHistory,costMatrix]},
adjustedForInflationCostTrajectory[notAdjusted,inflation]
],
CompilationTarget->"C",
CompilationOptions->{"InlineCompiledFunctions"->True,"ExpressionOptimization"->True},RuntimeOptions->"Speed"
];

costTrajectory[wh,𝒞,0.01]

Row[{Labeled[MatrixPlot[AntiDiagonalUpperTriangularToUpperTriangular[Module[{fasYears=5,colaRate=0.,employeeContributionRate=0.064},Table[PadRight[bothSides[salaryScale,fasYears,survival,survivalTableInitialAge,discountRate,inflation,txReplRate,True,True,colaRate,0.`,1000,killOffAge,{{0,1000000,employeeContributionRate}},employerContributionRatesList,serviceDomain,retirementDomain,{entryAge,0},multipleHazardTX[[entryAge-19]],discountRate,1,1][[-1]],56],{entryAge,20,75,1}]]
],FrameTicks->None],"normal cost matrix"],"×",Labeled[MatrixPlot[workforceTensor[[14]],FrameTicks->None],"workforce distribution\n
in year 13"]},Spacer[12]]

With[{m=AntiDiagonalUpperTriangularToUpperTriangular[Module[{fasYears=5,colaRate=0.,employeeContributionRate=0.064},Table[PadRight[bothSides[salaryScale,fasYears,survival,survivalTableInitialAge,discountRate,inflation,txReplRate,True,True,colaRate,0.`,1000,killOffAge,{{0,1000000,employeeContributionRate}},employerContributionRatesList,serviceDomain,retirementDomain,{entryAge,0},multipleHazardTX[[entryAge-19]],discountRate,1,1][[-1]],56],{entryAge,20,75,1}]]
]},ListAnimate[MapIndexed[{matrix,index}|->Row[{Labeled[MatrixPlot[m,FrameTicks->None],"normal cost matrix"],"×",Labeled[MatrixPlot[matrix,FrameTicks->None],"workforce distribution\n
in year "<>ToString[index[[1]]-1]],Total[Flatten[m*matrix]]},Spacer[12]],workforceTensor]
]
]

Manipulate[Module[{yearsToForecast=50,workforceHistory,costMatrix,annualNormalCost,annualSalary,result},
costMatrix=Table[PadRight[Thread[bothSides[salaryScale,fasYears,survival,survivalTableInitialAge,discountRate,inflation,txReplRate,True,True,colaRate,0.`,1000,killOffAge,{{0,1000000,employeeContributionRate}},employerContributionRatesList,serviceDomain,retirementDomain,{entryAge,0},multipleHazardTX[[entryAge-19]],discountRate,1,1][[{1,-1}]]],56,item[0.001,0]],{entryAge,20,75,1}];
workforceHistory=workforceEvolve[1-multipleHazardTX,workforceSnapshot,((Diagonal[workforceSnapshot])/(Total[Diagonal[workforceSnapshot]])),growthRate,yearsToForecast];
annualNormalCost=costTrajectory[workforceHistory,Map[Last,costMatrix,{2}],inflation];
annualSalary=costTrajectory[workforceHistory,Map[First,costMatrix,{2}],inflation];


result=annualNormalCost/annualSalary;
Labeled[ListPlot[Switch[plotWhat,"normal cost",annualNormalCost,"payroll",annualSalary,"normal cost percentage",result],
Joined->True,ImageSize->500,PlotMarkers->{Automatic,Tiny},Frame->True,Axes->False,FrameLabel->{"years into the future",plotWhat},GridLines->{Automatic,Automatic},ImagePadding->{{70,20},{60,10}},BaseStyle->{FontSize->14,FontFamily->"Swiss"}],
Style[plotWhat,{FontSize->14,FontFamily->"Swiss"}]]
],
{{plotWhat,"normal cost","plot what?"},{"normal cost","payroll","normal cost percentage"}},
Row[{"initial pension annuity is replacement rate multiplied by the average salary for the final ... ",Control@{{fasYears,5,""},Range[7],ControlType->SetterBar},"... years"}],
{{discountRate,0.05,"discount rate"},0.01,0.1,0.005,Appearance->"Labeled"},
{{inflation,0.025,"inflation rate"},0.,0.1,0.005,Appearance->"Labeled"},
{{colaRate,0,"COLA Rate"},0,0.1,0.005,Appearance->"Labeled"},
{{employeeContributionRate,0.064,"employee contribution rate"},0.0,0.08,0.001,Appearance->"Labeled"},
Delimiter,
{{growthRate,0,"rate of growth of workforce"},-0.02,0.04,0.001,Appearance->"Labeled"},
LabelStyle->{FontSize->14,FontFamily->"Swiss"}

]

Other materials for diagrams

Labeled[With[{entryAge=34},ListLinePlot[DeleteCases[Thread[List[Range[34,75],bothSides[salaryScale,fasYears,survival,survivalTableInitialAge,discountRate,inflation,txReplRate,True,True,0.,0.`,1000,killOffAge,{{0,1000000,0.064}},employerContributionRatesList,serviceDomain,retirementDomain,{entryAge,0},multipleHazardTX[[entryAge-19]],discountRate,1,1][[5]]]],{_,75.}],PlotMarkers->Automatic,Axes->False,Frame->True,FrameLabel->{"separation age","optimal retirement age"},ImageSize->720,BaseStyle->{FontSize->16,FontFamily->"Swiss"},PlotRange->{{34,75},{52,75}},GridLines->Automatic,Epilog->{Dashed,Gray,Opacity[0.6],Rectangle[{34,0},{39,100}],Opacity[1],Line[{{39,0},{39,100}}],Black,Text["refund\nzone",{35,70},{-1,0}]}]
],Style["Sample optimal retirement age for each separation age for a person with entry age of 34",{FontSize->16,FontFamily->"Swiss"}]]

Labeled[With[{entryAge=34},ListLinePlot[Thread[List[Range[34,75],bothSides[salaryScale,fasYears,survival,survivalTableInitialAge,discountRate,inflation,txReplRate,True,True,0.,0.`,1000,killOffAge,{{0,1000000,0.064}},employerContributionRatesList,serviceDomain,retirementDomain,{entryAge,0},multipleHazardTX[[entryAge-19]],discountRate,1,1][[-9]]]],PlotMarkers->Automatic,Axes->False,Frame->True,FrameLabel->{"separation age","present value\nat time of separation of pension benefit"},ImageSize->720,BaseStyle->{FontSize->16,FontFamily->"Swiss"},PlotRange->{{34,75},Automatic},GridLines->Automatic,Epilog->{Dashed,Gray,Opacity[0.6],Rectangle[{34,0},{39,100}],Opacity[1],Line[{{39,0},{39,100}}],Black,Text["refund\nzone",{35,70},{-1,0}]}]
],Style["Sample optimal present value at time of separation\nof pension benefit for a person with entry age of 34",{FontSize->16,FontFamily->"Swiss"}]]

Labeled[With[{entryAge=34},ListLinePlot[Thread[List[Range[34,75],bothSides[salaryScale,fasYears,survival,survivalTableInitialAge,discountRate,inflation,txReplRate,True,True,0.,0.`,1000,killOffAge,{{0,1000000,0.064}},employerContributionRatesList,serviceDomain,retirementDomain,{entryAge,0},multipleHazardTX[[entryAge-19]],discountRate,1,1][[-9]]]],PlotMarkers->Automatic,Axes->False,Frame->True,FrameLabel->{"separation age","present value\nat time of separation of pension benefit"},ImageSize->720,BaseStyle->{FontSize->16,FontFamily->"Swiss"},PlotRange->{{34,75},Automatic},GridLines->Automatic,Epilog->{Dashed,Gray,Opacity[0.6],Rectangle[{34,0},{39,100}],Opacity[1],Line[{{39,0},{39,100}}],Black,Text["refund\nzone",{35,70},{-1,0}]}]
],Style["Sample optimal present value at time of entry\nof gross pension benefit for a person with entry age of 34",{FontSize->16,FontFamily->"Swiss"}]]

NSPD=Expectation[valueOfPensionBenefits[separationTime],separationTime\[Distributed]separationDistribution]

Solve[NSPD==Expectation[TimeValue[Cashflow[Table[{t,augmentedNormalCostPercentage*salary[t]},{t,0,separationTime-1}]],r,0],
separationTime\[Distributed]separationDistribution],
augmentedNormalCostPercentage]

Solve[NSPD==Expectation[TimeValue[Cashflow[Table[{t,normalCostPercentage*salary[t]},{t,0,separationTime-1}]],r,0],
separationTime\[Distributed]separationDistribution]+Expectation[TimeValue[Cashflow[Table[{t,employeeContribution[salary[t]]*salary[t]},{t,0,separationTime-1}]],r,0],
separationTime\[Distributed]separationDistribution],
normalCostPercentage]