Replicating the AAA equtity generator

using CairoMakie
using ColorSchemes
using Distributions
using LabelledArrays
using Random

This notebook replicates the model and parameters for the real world equity generator described in this AAA 2005 reference paper.

Stochastic Log Volatility Model

Note that the @. and other broadcasting (. symbol) allows us to operate on multiple funds at once.

function v(v_prior,params,Zₜ) 
    (;σ_v, σ_m,σ_p,σ⃰,ϕ,τ) = params
    
    v_m = log.(σ_m)
    v_p = log.(σ_p)
    v⃰ = log.(σ⃰)

    # vol are the odd values in the random array
    ṽ =  @. min(v_p, (1 - ϕ) * v_prior + ϕ * log(τ) ) + σ_v * Zₜ[[1,3,5,7]]
    
    v = @. max(v_m, min(v⃰,ṽ))

    return v
end

function scenario(params,Z;months=1200)
    (;σ_v,σ_0, ρ,A,B,C) = params

    n_funds = size(params,2)
    
    #initilize/pre-allocate
    Zₜ = rand(Z)
    v_t = log.(σ_0)
    σ_t = zeros(n_funds)
    μ_t = zeros(n_funds)
    
    log_returns = map(1:months) do t
        Zₜ = rand!(Z,Zₜ)
        v_t .= v(v_t,params,Zₜ)

        σ_t .= exp.(v_t)

        @. μ_t =  A + B * σ_t + C * (σ_t)^2

        # return are the even values in the random array
        log_return = @. μ_t / 12 + σ_t / sqrt(12) * Zₜ[[2,4,6,8]]
    end

    # convert vector of vector to matrix
    reduce(hcat,log_returns)
end

scenario (generic function with 1 method)

Model Parameters

# use a labelled array for easy reference of the parameters 
params = @LArray [
    0.12515 0.14506 0.16341 0.20201     # τ
    0.35229 0.41676 0.3632 0.35277      # ϕ
    0.32645 0.32634 0.35789 0.34302     # σ_v
    -0.2488 -0.1572 -0.2756 -0.2843     # ρ
    0.055 0.055 0.055 0.055             # A
    0.56 0.466 0.67 0.715               # B
    -0.9 -0.9 -0.95 -1.0                # C
    0.1476 0.1688 0.2049 0.2496         # σ_0
    0.0305 0.0354 0.0403 0.0492         # σ_m
    0.3 0.3 0.4 0.55                    # σ_p
    0.7988 0.4519 0.9463 1.1387         # σ⃰
] ( 
    # define the regions each label refers to
    τ = (1,:),
    ϕ = (2,:),
    σ_v = (3,:),
    ρ = (4,:),
    A = (5,:),
    B = (6,:),
    C = (7,:),
    σ_0 = (8,:),
    σ_m = (9,:),
    σ_p = (10,:),
    σ⃰ = (11,:)
)

11×4 LArray{Float64, 2, Matrix{Float64}, (τ = (1, Colon()), ϕ = (2, Colon()), σ_v = (3, Colon()), ρ = (4, Colon()), A = (5, Colon()), B = (6, Colon()), C = (7, Colon()), σ_0 = (8, Colon()), σ_m = (9, Colon()), σ_p = (10, Colon()), σ⃰ = (11, Colon()))}:
   :τ => 0.12515  …    :τ => 0.20201
   :ϕ => 0.35229       :ϕ => 0.35277
 :σ_v => 0.32645     :σ_v => 0.34302
   :ρ => -0.2488       :ρ => -0.2843
   :A => 0.055         :A => 0.055
   :B => 0.56     …    :B => 0.715
   :C => -0.9          :C => -1.0
 :σ_0 => 0.1476      :σ_0 => 0.2496
 :σ_m => 0.0305      :σ_m => 0.0492
 :σ_p => 0.3         :σ_p => 0.55
   :σ⃰ => 0.7988   …    :σ⃰ => 1.1387

The Multivariate normal and covariance matrix

# 11 columns because it's got the bond returns in it
cov_matrix = [
    1.000   -0.249  0.318   -0.082  0.625   -0.169  0.309   -0.183  0.023   0.075   0.080;
    -0.249  1.000   -0.046  0.630   -0.123  0.829   -0.136  0.665   -0.120  0.192   0.393;
    0.318   -0.046  1.000   -0.157  0.259   -0.050  0.236   -0.074  -0.066  0.034   0.044;
    -0.082  0.630   -0.157  1.000   -0.063  0.515   -0.098  0.558   -0.105  0.130   0.234;
    0.625   -0.123  0.259   -0.063  1.000   -0.276  0.377   -0.180  0.034   0.028   0.054;
    -0.169  0.829   -0.050  0.515   -0.276  1.000   -0.142  0.649   -0.106  0.067   0.267;
    0.309   -0.136  0.236   -0.098  0.377   -0.142  1.000   -0.284  0.026   0.006   0.045;
    -0.183  0.665   -0.074  0.558   -0.180  0.649   -0.284  1.000   0.034   -0.091  -0.002;
    0.023   -0.120  -0.066  -0.105  0.034   -0.106  0.026   0.034   1.000   0.047   -0.028;
    0.075   0.192   0.034   0.130   0.028   0.067   0.006   -0.091  0.047   1.000   0.697;
    0.080   0.393   0.044   0.234   0.054   0.267   0.045   -0.002  -0.028  0.697   1.000;
]

    Z = MvNormal(
        zeros(11), #means for return and volatility
        cov_matrix # covariance matrix
        # full covariance matrix in AAA Excel workook on Parameters tab
    )

FullNormal(
dim: 11
μ: [0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0]
Σ: [1.0 -0.249 … 0.075 0.08; -0.249 1.0 … 0.192 0.393; … ; 0.075 0.192 … 1.0 0.697; 0.08 0.393 … 0.697 1.0]
)

Scenarios and validation

A single scenario

x = scenario(params,Z;months=1200)

4×1200 Matrix{Float64}:
 0.0752418  -0.000684304  0.0803535  …   0.000348079  -0.00459584
 0.0253968   0.0235449    0.0841353      0.0127518    -0.0323876
 0.0893632   0.0144473    0.0872761     -0.0191581     0.00114621
 0.138261    0.0343143    0.0945716     -0.0176436    -0.029718

Validation of summary statistics

The summary statistics expected (per paper Table 8):

μ ≈ [0.0060, 0.0062, 0.0063, 0.0065]
σ ≈ [0.0436, 0.0492, 0.0590, 0.0724]

These computed values match very closely:

# generate 1000 scenarios 
scens = [scenario(params,Z) for _ in 1:1000];

let
    # compute summary statistics
    μ = vec(mean(mean(x,dims=2) for x in scens))
    σ = vec(mean(std(x,dims=2) for x in scens))
    (;μ,σ)
end

(μ = [0.005996167539900355, 0.006165861874564148, 0.006254386906381312, 0.006366151507437672], σ = [0.0436367414351783, 0.04920555060632976, 0.05914308496454576, 0.0724382319579281])

Plotting some scenarios

let 
    f = Figure()
    n = 25
    colors = ColorSchemes.Johnson
    ax = Axis(f[1,1],yscale=log10,ylabel="index value",
        title="$n realizations of 4 correlated equity funds per AAA ESG")
    for s in scens[1:n]
        for i in 1:4
        lines!(ax,cumprod(exp.(s[i,:])), color=(colors[i],0.3),label="fund $i")
        end
    end
    axislegend(ax,unique=true,position=:lt)
    f
end

┌ Warning: Found `resolution` in the theme when creating a `Scene`. The `resolution` keyword for `Scene`s and `Figure`s has been deprecated. Use `Figure(; size = ...` or `Scene(; size = ...)` instead, which better reflects that this is a unitless size and not a pixel resolution. The key could also come from `set_theme!` calls or related theming functions.
└ @ Makie ~/.julia/packages/Makie/We6MY/src/scenes.jl:227