Introduction
I was having a conversation with a colleague about the lack of performance in actuarial models the other day. We started talking about the Academy's Interest Rate Generator (AIRG). Was it performant? My colleague argued that it calculates 10,000 scenarios in twenty minutes. That's over eight scenarios a second! He could run it multiple times in a day if he needed! I was not impressed. but it's hard to convince an actuary without some numbers backing it up.
History
The American Academy of Actuaries’ Economic Scenario Working Group (ESWG) initially released its interest rate scenario generator model in 1999 as a part of the C-3 Phase I project. The largest update since then was in December of 2008 where the parameters were updated, other models were considered, and calibration criteria for a company's own generator were designed. Since 2008 it has gone through minor tweaks and enhancements, but the model remains in Excel – specifically VBA.
I've included some links below to more information on the model and its updates. It was frustrating to research the AIRG's history as none of it is in one comprehensive location. Hopefully these links save you the time that I'll never get back.
- C3-Phase 1 Report - October 1999
- NAIC 2007 Q3 Meeting - pages 1446 - 1448 and pages 1456 - 1464
- AAA December 2008 Report
- AAA ESWG March 2007 Presentation
- AAA ESWG September 2006 Presentation
- AAA ESG Version 7 Release Notes
- AAA AIRG FAQ May 2019
Disclaimer
First, you should absolutely not use the code below in any form. This is for demonstration purposes only.
Second, I understand that Excel and VBA aren't built for performance. However, that hasn't stopped actuaries from building critical work in Excel. Waiting for a long-winded Excel model to run isn't adding value. Actuaries can and should do better.
The Model
The AIRG is a stochastic log volatility model. This means that the volatility of a stochastic process is also stochastic.
The model defines three stochastic processes:
- The natural log of the long term interest rate at time t
- The excess of the short term rate over the long term rate at time t
- the natural log of the monthly variance of stochastic process 1
where t is a monthly time step and where the stochastic process is a standard normal random variable.
The Excel model uses a hand-written random number generator and rewrites an inverse-normal function in order to determine the processes. Julia's base library contains many RNG algorithms to avoid having to write our own. For this example we'll use the Mersenne Twister algorithm that the AIRG VBA makes a comment about.
function generate_normals(months::Int)
rng = MersenneTwister()
randn(rng,(3, months))
end
It's easy enough to pass in a seed to the RNG by utilizing Julia's multiple dispatch.
function generate_normals(months::Int, seed::Int)
rng = MersenneTwister(seed)
randn(rng,(3, months))
end
The model takes the stochastic processes and converts them to correlated normals using the following formulas:
\begin{align} _1Z_t & = _1N_t \\ _2Z_t & = _1N_t * \rho_{1,2} + _2N_t * \sqrt{1-\rho_{1,2}^2} \\ _3Z_t & = _1N_t * \rho_{1,3} + \frac{_2N_t * (\rho_{2,3} - \rho_{1,2} * \rho_{1,3})}{\sqrt{1-\rho_{1,2}^2}} + _3N_t * \sqrt{1 - \frac{(\rho_{2,3} - \rho_{1,2} * \rho_{1,3})^2}{1-\rho_{1,2}^2} - \rho_{1,3}^2} \end{align}
Did you recognize the formulas from the Cholesky decomposition? That's not referenced anywhere. By knowing this, we're able to leverage some high performance linear algebra code to do the hard work for us and avoid having to code all of the formulas by hand.
function correlate_normals!(N::Array{Float64,2}, ρ::Matrix{Float64})
lmul!(cholesky!(ρ).L,N)
end
The correlation matrix used to get Cholesky decomposition in the Excel model takes a little searching and reconstructing, but you end up with:
$$
\begin{bmatrix}
1 & -0.19197 & 0 \\
-0.19197 & 1 & 0 \\
0 & 0 & 1
\end{bmatrix}
$$
We'll use this later when we start running things, but for now let's store it off in a variable.
p = [1.0 -0.19197 0.0; -0.19197 1.0 0.0; 0.0 0.0 1.0]
Once we have our correlated normals, we're able to get to the meat of the work. The model is mean reverting for all three processes. It uses $\tau$ for each mean reversion point and $\beta$ for the strength of the mean reversion. They also set some bounds on the long interest rate. The formulas are as follows:
\begin{align} _1i_t & = max(_1\lambda_L, min(_1\lambda_U, (1-\beta_1) * _1i_{t-1} + \beta_1 * \ln \tau_1 + \Psi * (\tau_2 - \alpha_{t-1}))) + _1\sigma_t * _1Z_t \\ \alpha_t & = (1-\beta_2) * \alpha_{t-1} + \beta_2 * \tau_2 + \phi * (_1i_{t-1} - \ln \tau_1) + \sigma_2 * _2Z_t * _1r_{t-1}^\theta \\ \nu_t & = (1-\beta_3) * \nu_{t-1} + \beta_3 * \ln \tau_3 + \sigma_3 * _3Z_t \\ & where \\ _1i_t & = \ln _1r_t \\ _1\lambda_L & = \ln _1r_{min} \\ _1\lambda_U & = \ln _1r_{max} \\ _2r_t & = \exp(_1i_t) - \alpha_t \\ _1\sigma_t & = \exp(\nu_t) \end{align}
A lot of parameters end up getting used for these. By creating a function, we're able to pass in the AIRG's defaults.
function generate_processes(Z::Array{Float64,2},
init_r_short,
init_r_long,
r_min=0.015,
r_max=0.18
Β_1=0.00509,
Β_2=0.02685,
Β_3=0.04001,
τ_1=0.0350,
τ_2=0.01,
τ_3=0.0287,
σ_1=0.0287,
σ_2=0.04148,
σ_3=0.11489,
Ψ=0.25164,
Φ=0.0002,
θ=1.0
)
ν = Vector{Float64}(undef, size(Z,2)+1)
α = similar(ν)
i = similar(ν)
for j = 1:size(ν,1)
if j == 1
ν[j] = log(σ_1)
α[j] = init_r_long - init_r_short
i[j] = log(init_r_long)
else
ν[j] = (1.0- Β_3) * ν[j-1] + Β_3 * log(τ_3) + σ_3 * Z[3,j-1]
i[j] = max(log(r_min),min(log(r_max),(1 - Β_1) * i[j-1] + Β_1 * log(τ_1) + Ψ * (τ_2 - α[j-1]))) + exp(ν[j]) * Z[1,j-1]
α[j] = (1.0 - Β_2) * α[j-1] + Β_2 * τ_2 + Φ * (i[j-1] - log(τ_1)) + σ_2 * Z[2,j-1] * exp(i[j-1]) ^ θ
end
end
return i, α
end
From here we calculate our fitted parameters:
\begin{align} b_0 &= _2r - b_1 * ns_1 \\ b_1 &= \frac{_2r - _1r}{ns_1-ns_{20}} \end{align}
where ns are the Nelson-Siegel interpolation factors for year 1 and year 20 respectively:
$$ ns = \frac{1-\exp(-0.4*t)}{0.4*t} $$
and in Julia:
function fit_parameters(i, α, short_year, long_year)
short_year_factor = nelson_siegel(short_year)
long_year_factor = nelson_siegel(long_year)
b_0 = similar(i)
b_1 = similar(i)
@. b_1 = ((exp(i) - α) - exp(i)) / (long_year_factor - short_year_factor)
@. b_0 = exp(i) - α - b_1 * short_year_factor
return b_0, b_1
end
and
function nelson_siegel(t::Float64)
(1.0 - exp(-0.4 * t))/(0.4*t)
end
We use these Nelson-Siegel factors for all ten points along the yield curve and for all time $t$.
function interpolate_curve(b_0, b_1, interpolation_factors)
curve = Array{Float64}(undef, size(b_0,1), size(interpolation_factors,1))
@. curve = b_0 + b_1 * interpolation_factors'
return curve
end
The AIRG doesn't stop after interpolation though. It perturbs the interpolated values to align time zero with the actual yield curve – grading off over twelve months.
function perturb_curve!(fitted_yc, init_yc, months)
perturb = zeros(size(fitted_yc,1))
perturb[1:months+1] = collect(1:-1.0/months:0)
@. fitted_yc = (1.0 - perturb) * fitted_yc + perturb * init_yc'
end
Composing our functions results in our end products:
function generate_yieldCurve(yieldCurve, ρ, forecast_months, perturb_months)
tc = collect(keys(yieldCurve))
yc = collect(values(yieldCurve))
interpolation_factors = similar(tc)
@. interpolation_factors = nelson_siegel(tc)
proc = generate_processes(correlate_normals!(generate_normals(forecast_months), ρ))
params = fit_parameters(proc[1], proc[2], 1, 20)
curve = interpolate_curve(params[1], params[2], interpolation_factors)
perturb_curve!(curve,yc,perturb_months)
return curve
end
where the $yieldCurve$ variable is an ordered dictionary from the DataStructures package.
We can define our starting yield curve by:
term_structure = [0.25,0.5,1.0,2.0,3.0,5.0,7.0,10.0,20.0,30.0]
initYC = [0.0245,0.0256,0.0263,0.0248,0.0246,0.0251,0.0259,0.2069,0.0287,0.0302]
tc = OrderedDict(term_structure .=> initYC)
By using the BenchmarkTools package, we can get an idea of what kind of performance to expect and any potential bottlenecks:
@benchmark generate_yieldCurve(tc,p,360,12)
This macro samples the function 10,000 times and gives us some metrics:
BenchmarkTools.Trial:
memory estimate: 75.47 KiB
allocs estimate: 31
--------------
minimum time: 84.600 μs (0.00% GC)
median time: 99.900 μs (0.00% GC)
mean time: 103.411 μs (0.00% GC)
maximum time: 322.100 μs (0.00% GC)
--------------
samples: 10000
evals/sample: 1
I'm certain there are potential enhancements in my code as I'm not a Julia expert by any means, but I'll take 100 microseconds with no garbage collection for a scenario as a win.
Computing more than one scenario is as simple as writing a loop:
function generate_scenarios(yieldCurve, ρ, forecast_months, perturb_months, scenarios)
for i = 1:scenarios
generate_yieldCurve(tc,p,forecast_months,perturb_months)
end
end
and if we want to take advantage of the fact that each scenario is independent and can be run on separate threads, we can parallelize the scenario generation with the @threads macro:
function parallel_generate_scenarios(yieldCurve, ρ, forecast_months, perturb_months, scenarios)
Threads.@threads for i = 1:scenarios
generate_yieldCurve(tc,p,forecast_months,perturb_months)
end
end
Comparison
I tested both the Excel model and the Julia code on my i7 laptop.
- The Excel model generated 10,000 scenarios in 8 minutes and 26 seconds.
- The Julia code generated 10,000 scenarios in 1.157 seconds.
- The parallel Julia code generated 10,000 scenarios in 346.098 milliseconds.
The parallel Julia code is a 1,462x performance improvement over the VBA code.
Conclusion
This exercise was more frustrating than it should have been due to the AIRG's documentation (or lack there of). All of the calculations are done in VBA with the Excel sheet only providing the calculations post run. I ran into multiple occasions where formulas or calculations were either deprecated or not used, and the documentation failed to be updated to reflect that.
The Julia code in this article is by no means perfectly written, optimized, or feature complete with the AIRG Excel model. Even still, the performance improvements hopefully demonstrate that their are many opportunities within the actuarial profession to become more efficient and add more value.
What are your thoughts on performance and efficiency in actuarial work? Do you agree with my assessment? What are the low hanging fruit in the industry? Email me any comments or questions you may have at hp@houstonp.com.