4.1 Additive and multiplicative time series

In Derman & Kani (Reference Derman and Kani1994), Dupire (Reference Dupire1994), and Rubinstein (Reference Rubinstein1994), it is postulated that the volatility of assets is not only time-inhomogeneous, but it is also a function of the prices of the assets themselves or, by extension, the recent changes in the asset prices – this is inline with modern probability default models that agree that the momentum of the underlying factors is nontrivial in their predictive capabilities. This follows logically as high prices or recent large increases in asset prices often foreshadow periods of heightened uncertainty and, therefore, volatility. These assumptions are not without criticism with Dumas et al. (Reference Dumas, Fleming and Whaley1998), performing an extensive empirical study that concludes that these models are indistinguishable from smoothed versions of the implied volatility model proposed in Black & Scholes (Reference Black and Scholes1973).

With this context, it is natural to also explore multiplicative models where seasonality and trend components (and further less-easily named or interpretable components isolable using SSA and SSD within the CovRegpy package and EMD in the AdvEMDpy package) are multiplied onto some underlying trend. The framework in this paper is also well suited for multiplicative models, provided that one uses an appropriate transform such as the natural logarithm as demonstrated in Equation (6):

(6) \begin{equation} \begin{aligned} x_t &= x_{t, \text{trend}} \times x_{t, \text{seasonal}} \times x_{t, \text{irregular}}\\ \text{log}(x_t) &= \text{log}(x_{t, \text{trend}}) + \text{log}(x_{t, \text{seasonal}}) + \text{log}(x_{t, \text{irregular}}). \end{aligned} \end{equation}

4.2 Empirical mode decomposition

Based on work done in Huang et al. (Reference Huang, Shen, Long, Wu, Shih, Zheng, Yen, Tung and Liu1998), Huang (Reference Huang1999), Huang et al. (Reference Huang, Shen and Long1999), van Jaarsveldt et al. (Reference van Jaarsveldt, Ames, Peters and Chantler2023a), and van Jaarsveldt et al. (Reference van Jaarsveldt, Peters, Ames and Chantler2023b), it can be stated with some confidence that EMD and the various algorithmic variations summarized and presented in van Jaarsveldt et al. (Reference van Jaarsveldt, Ames, Peters and Chantler2023a) present a class of robust trend estimation and decomposition algorithm that are state-of-the-art when working with general stochastic processes including potential nonstationarity and nonlinearity. The algorithm, Algorithm 1, presented in pseudocode in Appendix A, of Supplement to: Package CovRegpy: Regularized Covariance Regression and Forecasting in Python, was first presented in van Jaarsveldt et al. (Reference van Jaarsveldt, Ames, Peters and Chantler2023a). The two conditions that each component isolated using EMD, named intrinsic mode functions (IMFs), must satisfy are listed in C1 and C2. An additional, optional (and highly recommended) condition, C3, is referred to as a stopping criterion and is checked after the initial two conditions are checked. This prevents over-sifting and the propagation of errors through successive IMF components with the conditions listed below:

1. C1 $ \text{abs}\bigg (\Big |\Big \{\frac{d\gamma _k(t)}{dt} = 0 \,{:}\, t \in (0, T)\Big \}\Big | - \Big |\Big \{\gamma _k(t) = 0 \,{:}\, t \in (0, T)\Big \}\Big |\bigg ) \leq 1,$

2. C2 $ \sum _t\text{abs}\big (\tilde{\gamma }^{\mu }_k(t) \big ) \leq \epsilon _2, \text{ and}$

3. C3 $ SD_{(p,q)} = \sum _t \Bigg [\dfrac{\big |(h_{(p,q-1)}(t) - h_{(p,q)}(t))\big |^2}{h^2_{(p,q-1)}(t)}\Bigg ] \lt \epsilon _3.$

In C1, $ \gamma _k(t)$ is the k $^{th}$ IMF component and with $ |\{\cdot \}|$ being the cardinality set. In C2, $ \tilde{\gamma }^{\mu }_k(t)$ is the mean of the k $^{th}$ IMF and in C3, $ h_{(p,q)}(t)$ is the q $^{th}$ iteration of the algorithm for the p $^{th}$ IMF. The final result of the sifting presented in Algorithm 1 is

(7) \begin{equation} x_{\text{IMF}}(t) = \sum _{k=1}^{K}\gamma _k(t) + r_K(t). \end{equation}

As mentioned, EMD refers exclusively to the sifting procedure, whereas the Hilbert–Huang transform (HHT) refers to both EMD (the sifting procedure) followed by the Hilbert transform. The Hilbert transform of $\textrm{kth}$ IMF can be seen in Equation (8):

(8) \begin{equation} \check{\gamma }_k(t) = \mathcal{HT}[\gamma _k(t)] = \frac{1}{\pi }PV\int _{-\infty }^{\infty }\dfrac{\gamma _k(t^*)}{t - t^*}dt^*, \end{equation}

with $PV$ being the Cauchy principle value integral. This is an improper integral that assigns a value to an improper integral where there exists a discontinuity. The Hilbert transform in Equation (8) can be expanded to Equation (9):

(9) \begin{equation} \check{\gamma }_k(t) = \frac{1}{\pi }\lim _{\epsilon \rightarrow{0}^+}\bigg [\int _{t - \frac{1}{\epsilon }}^{t - \epsilon }\dfrac{\gamma _k(t^*)}{t - t^*}dt^* + \int _{t^* + \epsilon }^{t + \frac{1}{\epsilon }}\dfrac{\gamma _k(t^*)}{t - t^*}dt^*\bigg ]. \end{equation}

With the Hilbert transform of the time series defined as in Equation (9), the analytical signal can then be defined as in Equation (10):

(10) \begin{equation} \gamma ^a_k(t) = \gamma _k(t) + i\check{\gamma }_k(t) = a(t)e^{i\theta (t)}= a(t)e^{i\int \omega (t)dt}, \end{equation}

with the instantaneous amplitude being defined as in Equation (11):

(11) \begin{equation} a(t) = \sqrt{\gamma _k(t)^2 + \check{\gamma }_k(t)^2}, \end{equation}

and the phase being:

(12) \begin{equation} \theta (t) = \text{tan}^{-1}\left (\frac{\check{\gamma }_k(t)}{\gamma _k(t)}\right ). \end{equation}

The instantaneous frequency can then be defined as in Equation (13):

(13) \begin{equation} \omega (t) = \frac{d\theta (t)}{dt}. \end{equation}

4.3 X11

X11 was originally proposed in Shiskin et al. (Reference Shiskin, Young and Musgrave1967). This algorithm is the focus of CovRegpy_X11.py. A disadvantage of X11 when compared against EMD is that it is less flexible in that it can only extract three structures namely the trend-cycle component, seasonal component, and the irregular component. Each of these components is individually less flexible than an IMF from EMD. The irregular component would, in most real-world scenarios, contain additional meaningful structural information that is regarded as noise in X11. Many shortcomings of the method are discussed in Sutcliffe (Reference Sutcliffe1993) and Doherty (Reference Doherty2001).

The Henderson symmetric weights proposed in Henderson (Reference Henderson1916), Whittaker (Reference Whittaker1922), and Henderson (Reference Henderson1924) are used in the X11 algorithm. The edges of the time series present problems for the symmetric weighting procedure. Two common solutions are proposed. As in Dagum (Reference Dagum1980), (Reference Dagum1988), (Reference Dagum1996), and Findley et al. (Reference Findley, Monsell, Bell, Otto and Chen1998), the autoregressive integrated moving average (ARIMA) model is used to explicitly forecast the edge of the time series so that the symmetric weights can still be used. An alternative to this without explicitly forecasting the time series involves calculating a set of asymmetric weights. Classical asymmetric weights are calculated in Musgrave (Reference Musgrave1964b) and (Reference Musgrave1964a) with the Reproducing Kernel Hilbert Space derivation of these weights being done Bianconcini (Reference Bianconcini2006), Dagum & Bianconcini (Reference Dagum and Bianconcini2006), and (Reference Dagum and Bianconcini2008). X11 decomposes the times series as:

(14) \begin{equation} x(t) = T(t) + S(t) + \epsilon (t), \end{equation}

with $ T(t)$ , $ S(t)$ , and $ \epsilon (t)$ being the trend-cycle component, seasonal component, and irregular component, respectively. The symmetric moving average of order $ k$ , $ \text{ma}_k(\cdot )$ , is defined as:

(15) \begin{equation} \text{ma}_k(x(t_h)) = \frac{1}{2(k-1)}x\Big (t_{h-\frac{(k-1)}{2}}\Big ) + \mathop{\sum }\limits _{a=-\frac{(k-2)}{2}}^{\frac{(k-2)}{2}}\frac{1}{k-1}x(t_{h+a}) + \frac{1}{2(k-1)}x\Big (t_{h+\frac{(k-1)}{2}}\Big ), \end{equation}

with $ x(t_h)$ being the $ h^{th}$ time series component. The symmetric seasonal moving average of order $ m\times{n}$ , $ S_{\{{m\times{n}}, c\}}(\cdot )$ , is defined as:

(16) \begin{equation} \begin{aligned} S_{\{{m\times{n}}, c\}}(x(t_h)) &= \frac{1}{mn}x\Big (t_{h-\frac{c((n-1)+(m-1))}{2}}\Big ) + \frac{2}{mn}\mathop{\sum }\limits _{b=-\frac{c((n-1)+(m-1)-2)}{2}}^{-\frac{c(n-1)}{2}}x(t_{h+b})\\ &+ \mathop{\sum }\limits _{b=-\frac{c(n-3)}{2}}^{\frac{c(n-3)}{2}}\frac{1}{n}x(t_{h+b}) + \frac{2}{mn}\mathop{\sum }\limits ^{\frac{c((n-1)+(m-1)-2)}{2}}_{b=\frac{c(n-1)}{2}}x(t_{h+b})\\ &+ \frac{1}{mn}x\Big (t_{h+\frac{c((n-1)+(m-1))}{2}}\Big ), \end{aligned} \end{equation}

with $ c$ being a parameter that depends on both the sampling rate of the time series and the calculated seasonality of the time series. An example is that if the time series is sampled monthly and the seasonality is expected to be quarterly, then $ c = 3$ . The Henderson moving average of order $ l$ , $ \text{Hma}_l(\cdot )$ , is as in Henderson (Reference Henderson1916), Whittaker (Reference Whittaker1922), Henderson (Reference Henderson1916), Bianconcini (Reference Bianconcini2006), Dagum & Bianconcini (Reference Dagum and Bianconcini2006), and Dagum & Bianconcini (Reference Dagum and Bianconcini2008). An example of the application of the code can be seen below:

trend, seasonal, irregular = \

CovRegpy_X11(time, time_series, seasonality=‘annual’,

seasonal_factor=‘3x3’, trend_window_width_1 = 13,

trend_window_width_2 = 13, trend_window_width_3 = 13)

Fig. 3 demonstrates the result of applying the X11 method in CovRegpy_X11 to the additive synthetic time series, Fig. 2, described in Section 4 using the parameters as shown in the code above. The trend-cycle component, seasonal component, and random error component are separated from one another, but the noise persists in the system owing to the lack of formal frequency bandwidth or distribution approaches to the decomposition as will be formally demonstrated in Section 4.5. The three components visible in Fig. 3 result from applying the above code snippet to the time series presented in Equation (2) and plotting each component with its corresponding underlying target time series component.

Figure 3 X11 time series decomposition of synthetic additive time series into trend-cycle component, seasonal component, and random error or noise component.

4.4 Singular spectrum analysis

SSA, as presented in Hassani (Reference Hassani2007), can be separated into four stages, namely embedding, singular value decomposition (SVD), grouping, and diagonal averaging which will be presented in Sections 4.4.1, 4.4.2, 4.4.3, and 4.4.4, respectively. In Section 9.1.1 of Supplement to: Package CovRegpy: Regularized Covariance Regression and Forecasting in Python, D-SSA is introduced and in Section 9.1.2 of “Supplement to: Package CovRegpy: Regularized Covariance Regression and Forecasting in Python,” KS-SSA is introduced. D-SSA and KS-SSA have been developed during this research to modify SSA to develop a more robust decomposition technique. D-SSA seeks to construct a decomposing version of SSA, whereas KS-SSA was developed to optimize and automate the isolation of a trend among white noise – the assumption underpinning this technique is that the errors have a Gaussian distribution.

4.4.1 Embedding in SSA

The first step involves converting the univariate time series into a multivariate time series by embedding lagged increments of the original univariate time series to create matrix $ \boldsymbol{X}$ such that:

(17) \begin{equation} \boldsymbol{X} = [\boldsymbol{X}_1(t), \dots, \boldsymbol{X}_K(t)], \end{equation}

with $ \boldsymbol{X}_j(t) = [x(t_j), \dots, x(t_{j + L - 1})]^T$ and $ K = T - L + 1$ . $L$ is the embedding dimension that determines the structures being isolated and becomes increasingly relevant in the algorithmic variations and methodological extension of SSA, namely SSD (discussed in the next section), for purposes such as removing the initial trend structure.

4.4.2 Singular value decomposition in SSA

This is followed by the SVD of the matrix $ \boldsymbol{X}$ resulting in $ \boldsymbol{\lambda } = (\lambda _1, \dots, \lambda _L)$ and $ \textbf{U} = [\textbf{U}_1, \dots, \textbf{U}_L]$ being the vector of eigenvalues in descending order ( $ \lambda _1 \lt \dots \lt \lambda _L$ ) and their corresponding eigenvectors, respectively. With $ \textbf{V}_j ={\boldsymbol{X}^T}\textbf{U}_j/\sqrt{\lambda _j}$ , $ \boldsymbol{X}$ is decomposed as:

(18) \begin{equation} \boldsymbol{X} = \sum ^L_{j=1}{\boldsymbol{X}_j}, \end{equation}

where $ \boldsymbol{X}_j = \sqrt{\lambda }_j{\textbf{U}_j}{\textbf{V}^T_j}$ . Equation (18) demonstrates the decomposition of the embedding matrix in Equation (17). In the absence of rounding errors, Equation (18) should precisely reproduce the time series after a reversal of the embedding step.

4.4.3 Grouping in SSA

With $ \boldsymbol{I}$ being a subset of indices such that $ \boldsymbol{I} = \{i_1, \dots, i_p\} \subset \{1, \dots, L\} = \boldsymbol{L}$ , the trend is estimated using a user-selected subset of $ \boldsymbol{L}$ , but in practice to estimate the trend, owing to the exponential decay of the eigenvalue, only the first few indices are used:

(19) \begin{equation} \boldsymbol{X}_{\boldsymbol{I}} = \sum _{\boldsymbol{I}} \boldsymbol{X}_j. \end{equation}

Theoretically, any subset in $ \mathcal{P}(\boldsymbol{L})$ can be used to estimate the trend, but as the eigenvalues monotonically decrease, the amount of variation, and therefore information, in each subsequent $ \boldsymbol{X}_j$ decreases exponentially. It is this step that is adjusted in Section 9.1.1 of “Supplement to: Package CovRegpy: Regularized Covariance Regression and Forecasting in Python” to produce the D-SSA by keeping the components separate before performing the following step on individual components before grouping them based on the user’s requirements such as a pure trend estimate or an evaluation of different cyclical components.

4.4.4 Diagonal averaging in SSA

This final step is, in summary, a reversal of the embedding step. By appropriately defining a vector, such that the lagged embedding of Equation (17) is reversed so that an appropriate averaging can take place. With this vector, $ X_{I}[i]$ , defined as in Equation (20):

(20) \begin{equation} X_{I}[i] = \sum _{r+c=i}\boldsymbol{X}_{I}[r,c], \end{equation}

with the trend estimation then proceeding by using an appropriate averaging vector to reverse the lagged embedding:

(21) \begin{equation} x_{SSA}(t) = \sum _{i=1}^{p-1}\frac{(p+1)-i}{p}X_{I}[i] + \sum _{i=p}^{L-p}\frac{1}{p}X_{I}[i] + \sum _{i=L-p+1}^{L}\frac{(L-i)}{p}X_{I}[i]. \end{equation}

In the original SSA, as proposed in Hassani (Reference Hassani2007), the only input variables used in the code above are time_series, L, and est. The time series (time_series) is required as well as the embedding dimension (L) and the grouping factor (est). One variation on the original algorithm can be noted as the second output of the algorithm which is time_series_decomp. Rather than automatically group the structures, as in time_series_est, time_series_decomp stores each of the structures sequentially for the user’s convenience.

time_series_est, time_series_decomp = \

CovRegpy_ssa(time_series, L, est = 3, plot=False, KS_test=False,

plot_KS_test=False, KS_scale_limit = 1,

figure_plot=False)

An example of the D-SSA output (time_series_decomp) is plotted in Fig. 4. One is advantaged in this synthetic example in that one has foreknowledge of the underlying structures, but this would not necessarily hold when applied to real-world examples. The various components extracted can be seen against their corresponding underlying components in Fig. 4. The additional inputs in the above code will be discussed in the Supplement to: Package CovRegpy: Regularized Covariance Regression and Forecasting in Python.

Figure 4 SSA and D-SSA trend estimate and component isolation compared against the corresponding underlying structures.

4.5 Singular spectrum decomposition

In this section SSD, originally introduced in Bonizzi et al. (Reference Bonizzi, Karel, Meste and Peeters2014), is discussed and several algorithmic variations are also introduced. This technique is available to the user in CovRegpy_SSD.py as well as the algorithmic variations. For a detailed description of the modifications made to SSA to create SSD, see Bonizzi et al. (Reference Bonizzi, Karel, Meste and Peeters2014). In summation, the major adjustments are the creation and automation of the extraction of a trend described in Section 4.5.1 and the formal extraction of components based on defined narrow bandwidths as described in Section 4.5.2.

4.5.2 Downsampling in singular spectrum decomposition

This process can be summarized as calculating a frequency band which will be used to downsample the possibly detrended time series. By constructing three Gaussian functions such that:

(22) \begin{equation} \gamma (f, \theta ) = \sum _{i=1}^3{A_i}e^{-\frac{(f-\mu _i)^2}{2\sigma _i^2}}, \end{equation}

with the $ \mu$ parameters, and therefore, the centers of the Gaussian functions to be fitted to the power spectral density (PSD), being fixed and defined as:

(23) \begin{equation} \mu _1 = f_{\text{max}},\text{ }\mu _2 = f_2,\text{ }\mu _3 = \frac{f_{\text{max}} + f_2}{2}, \end{equation}

where $ f_{\text{max}}$ is the frequency at which the maximum PSD mode is located and $ f_2$ is the frequency at which the second-highest PSD mode is located, with the other parameters being initialized as:

(24) \begin{equation} \begin{array}{lcl} A^{(0)}_1 = \frac{1}{2}PSD(f_{\text{max}}), &\text{ }&\sigma ^{(0)}_1 = f\,{:}\,PSD(f) = \frac{2}{3}PSD(f_{\text{max}}),\\ A^{(0)}_2 = \frac{1}{2}PSD(f_2), &\text{ }&\sigma ^{(0)}_2 = f\,{:}\,PSD(f) = \frac{2}{3}PSD(f_2),\\ A^{(0)}_3 = \frac{1}{4}PSD(f_3), &\text{ }&\sigma ^{(0)}_3 = 4|f_{\text{max}} - f_2|, \end{array} \end{equation}

where $ \sigma ^{(0)}_1 = f\,{:}\,PSD(f) = \frac{2}{3}PSD(f_{\text{max}})$ denotes the frequency at which the PSD is equal to two-thirds of the PSD of the maximum PSD mode nearest to the maximum PSD – this allows for probable asymmetric peaks in the PSD.

L1 Singular Spectrum Decomposition:

With $ \mu _1$ , $ \mu _2$ , and $ \mu _3$ fixed, Equation (22) can be fitted using LASSO regression as in Equation (25). In practice, L1 SSD can result in a perpetual loop if not managed correctly and should be used with caution.

(25) \begin{equation} \text{min}_{\theta } ||\gamma (f, \theta ) - \text{PSD}(f)||_1. \end{equation}

L2 Singular Spectrum Decomposition:

With $ \mu _1$ , $ \mu _2$ , and $ \mu _3$ fixed, Equation (22) can be fitted can be fitted using ridge regression as in Equation (26). This is the preferred method as LASSO regression promotes scarcity and in the event of the $ \theta$ vector being optimized to be identically zeros, without a counting stopping criterion, would result in an infinite loop:

(26) \begin{equation} \text{min}_{\theta } ||\gamma (f, \theta ) - \text{PSD}(f)||^2_2. \end{equation}

Once Equation (22) has been fitted to the PSD, the optimal frequency band is calculated as $ (f_{\text{max}} - \delta{f}, f_{\text{max}} + \delta{f})$ where $ \delta{f}$ is calculated as in Equation (27):

(27) \begin{equation} \delta{f} = 2.5\sigma _1^{opt}. \end{equation}

This frequency band is demonstrated in Fig. 6. Figures, such as Figs. 5 and 6, can be optionally output when running the code below if one has plot=True which outputs every incremental PSD fitting process.

4.5.3 Stopping criterion in singular spectrum decomposition

To prevent over-sifting (a word borrowed from EMD, but accurately summarizes the concept), one should introduce a stopping criterion that stops the algorithm after a certain number of components have been isolated. The SSD algorithm results in the following decomposition:

(28) \begin{equation} x_{SSD}(t) = \sum _{i=1}^{M}\tilde{g}_i(t) + v_{M+1}(t), \end{equation}

with $ \tilde{g}_i(t)$ being the $\textrm{ith}$ component and $ v_{M+1}(t)$ being the residual. This stopping criterion stops the algorithm when a certain percentage of the PSD of the time series has been accounted for in the formal components extracted. The algorithm stops when the normalized mean squared error (NMSE), calculated as follows:

(29) \begin{equation} NMSE_i = \sum ^{N}_{t=0}\frac{v_{(i+1)}^2(t)}{x^2(t)}, \end{equation}

drops below a certain percentage. A value of $ \alpha =0.01$ is recommended in Bonizzi et al. (Reference Bonizzi, Karel, Meste and Peeters2014) but is adjustable in this package. The time series (time_series) is the sole input required for the algorithm with the other additional inputs required being initial_trend_ratio, nmse_threshold, plot, and debug.

ssd_decomp = CovRegpy_ssd(time_series, initial_trend_ratio = 3,

nmse_threshold = 0.01, plot=False, debug=False)

The nmse_threshold value has been discussed in this section with initial_trend_ratio being discussed in Section 9.2.3 of Supplement to: Package CovRegpy: Regularized Covariance Regression and Forecasting in Python. The plot=True input variable results in numerous incremental plots such as Figs. 5 and 6. If debug=True, then each incremental value of the calculated $ NMSE_i$ will be printed for debugging purposes. The components resulting from SSD being applied to the synthetic time series described in Fig. 2 can be seen in Fig. 7 for comparison with Figs. 3 and 4.

Figure 5 Initialization of Gaussian functions, Equation (22), to be fitted to the PSD of the additive synthetic example presented in Fig. 2 for downsampling in SSD.

Figure 6 Fitting of Gaussian functions, Equation (22), to the PSD of the additive synthetic example presented in Fig. 2 for downsampling in SSD.

Figure 7 SSD of example time series using adjustable initial trend isolation window with the ubiquitous problem that also permeates EMD analysis known as the edge effect or frequency leakage between trend-cycle component and random error or irregular component.

4.6 Implicit factor context with regularized covariance regression

As detailed in Fig. 1, the three implicit factor extraction algorithms described in this section isolate components that may then be used as independent covariates ( $\mathbf{X}$ ) in the RCR framework described in the following section. The implicit factor algorithms are not a necessary step, but it is advised that if covariance forecasting is intended that the components used in Algorithm 1 be smoothed or subjected to some form of preprocessing to ensure the forecasted covariance evolves smoothly over time as in Figs. 8, 9, 10, and 11.

5.1 Early lung disease detection and insurance benefits

In Pyenson et al. (Reference Pyenson, Sander, Jiang, Kahn and Mulshine2012), the short-term costs and long-term benefits of screening for lung cancer in high-risk populations (aged 50–64 years in the USA) is demonstrated. The further reduced cost of the FEV screening (as an indicator for cancer or further, more-expensive, testing) of at-risk populations in the UK and elsewhere would greatly reduce the strain on national health institutions. The FEV measures examined herein, the research of Cho & Stout-Delgado (Reference Cho and Stout-Delgado2020) which notes a distinctive linear decrease in lung capacity with age, and our research which puts confidence bounds on the lung capacity as a function of age and size would all assist in cheap, early, and accurate measurement of reduced lung capacity that would reduce the future cost of cancer treatment. The research of Pyenson et al. (Reference Pyenson, Sander, Jiang, Kahn and Mulshine2012) strongly indicates that cheap and early screening for lung cancer would greatly decrease the future strain of healthcare costs on insurers.

Fitch et al. (Reference Fitch, Iwasaki, Pyenson, Plauschinat and Zhang2011) assesses the prevalence of COPD (in patients without coexisiting asthma) in a large group of patients (44,366) with a focus on determining a link between the severity of the condition and the level of adherence to prescribed medication. This work assists in determining the distributional severity of COPD (grouped into mild, moderate, severe, and very severe) among those making claims which allows insurers to make reasonable adjustments to their health plans based on this distribution. In Section 5.2, we make clear the diagnostic value of the FEV measurement for further tests or underlying lung degeneration which would also assist in grouping insurance clients into probabilistic groups based on likely severity of condition and mortality rates applications to policy pricing.

It is noted in Pyenson et al. (Reference Pyenson, Henschke, Yankelevitz, Yip and Dec2014) that lung cancer is, by a wide margin, the greatest cause of cancer deaths worldwide. Low-dose computed tomography (LDCT) is used to confirm the presence of lung cancer, but these tests are expensive and only conducted on populations deemed to be at risk based solely on smoking history and age – these two factors are, admittedly, strongly correlated with lung cancer instances. FEV and other inexpensive lung function tests are effective earlier indicators of reduced lung function for insurers to plan and allocate funds for future LDCT tests and potential treatments. Measures of lung capacity relative to age and size of an individual (which can be further stratified by gender) among a broader grouping of at-risk individuals such as younger populations in polluted cities or with vocations where exposure to hazardous airborne particulates can assist in early detection and further save insurers future health expenditures. The FEV (modeled herein) and other lung capacity tests (such as FEV per second compared to total FEV) are shown in Section 5.2 to be diagnostically valuable indicators of diseased lungs.

Dash & Grimshaw (Reference Dash and Grimshaw1993) assess the costs associated with positive cancer diagnoses from a different perspective: Dread Disease contracts. Dread Disease contracts are insurance policies that have become increasingly popular in the UK (Dash & Grimshaw, Reference Dash and Grimshaw1993) and payout predetermined lump sums upon the realization of certain mortal contingencies. As stated above, lung cancer is the largest cause of cancer deaths and any indicators of potential future lung deficiencies (such as FEV discussed herein) that could lead to cancer could significantly decrease the future cost to the underwriters of these Dread Disease contracts. Changes in marketing policies in the UK have lead to an increased interest in these contracts.

Efird et al. (Reference Efird, Landrine, Shiue, O’Neal, Podder, Rosenman and Biswas2014) investigates the stage of lung cancer diagnosis as it relates to racial groups and health insurance policies. It is noted that survival rates among lung cancer diagnosees are intimately linked with the stage at which it is first detected. It is noted that the type of insurance policy did not relate to the stage of the cancer diagnosis. Cheap and diagnostically valuable information such as early indicators of decreased lung capacity (FEV discussed herein) could assist in ameliorating this notably significant disparity between racial groups which would assist with the costs to insurers as a result of late diagnosis and expensive late-stage treatment compared to comparatively cheaper earlier stage treatments.

Yang et al. (Reference Yang, Beymer and Suen2019) notes the increasing life expectancies of those living with HIV and AIDS as a result of the provision of antiretrovirals. This decreased HIV-related mortality has had the unintended consequence of increased exposure to chronic illnesses, particularly lung diseases, which are often exacerbated by HIV, AIDS, or the associated antiretrovirals. It is noted in Yang et al. (Reference Yang, Beymer and Suen2019) that the prescence of HIV or AIDs in an individual’s system results symptoms that mimic an increased rate of biological aging which, when taking into account with this research and the known linear degradation in lung capacity with age (Cho & Stout-Delgado, Reference Cho and Stout-Delgado2020), can be used in the early detection of lung functional decay for the provision, by large private insurers (Yang et al., Reference Yang, Beymer and Suen2019), of funding for the relevant treatment.

5.2 Forced expiratory volume and mortality

As early as 1846 in Hutchinson (Reference Hutchinson1846) (which cites far older work), the measure of the lung capacity as a proxy for the normal function of the lungs was established as a viable diagnostic tool for determining the presence of lung disease. Specifically, this work focused on the spirometer which measures the rate of air expulsion from the lung as well as the total lung expulsion and uses these measurements (FEV per second and total FEV) against normally distributed values based on an individuals height (as a proxy of size and, therefore, expected lung capacity), gender, and age when making probabilistic predictions (which would be greatly assisted by our research) about an individuals health. These measurements and an individuals expected health incomes would greatly assist in calculating competitive insurance premiums.

Casanova et al. (Reference Casanova, Cote, de Torres, Aguirre-Jaime, Marin, Pinto-Plata and Celli2005) provides evidence that supports the correlative link between decreased FEV and the increased likelihood of death related to COPD which is caused by a number of serious medical conditions such as emphysema and chronic bronchitis. It is noted in this study, Casanova et al. (Reference Casanova, Cote, de Torres, Aguirre-Jaime, Marin, Pinto-Plata and Celli2005), that age is also a contributory factor in the decrease of lung capacity measured by the proxy of lung capacity, FEV. The study of the covariance between FEV and age in this work should be used in conjunction with work such as Dyer (Reference Dyer2012) which studies the decreases in lung capacity with age in the absence of confounding medical causes to better predict mortality outcomes.

Further, Dykstra et al. (Reference Dykstra, Scanlon, Kester, Beck and Enright1999) uses the FEV (among numerous other lung capacity measures) as a contributory diagnostic tool when determining whether a patient with a ratio of FEV per second to total FEV of less than 70% has asthma or COPD. By noting the stabilizing trend of FEV versus age in Fig. 8 and the work in Cho & Stout-Delgado (Reference Cho and Stout-Delgado2020) that supports a known approximately linear level of lung deterioration with the age beyond 35 years (approximately 20 mL/year), one can use deviations from the this trend as a contributory diagnostic factor. Cho & Stout-Delgado (Reference Cho and Stout-Delgado2020) and our work can also assist in the construction of modified mortality tables for individuals whose lung capacity depreciation rate is above those forecasted and observed.

In the study conducted in French et al. (Reference French, Balfe, Mirocha, Falk and Mosenifar2015), both age and lowered FEV (and the further measure of the ratio of inspiratory capacity to total lung capacity once a lowered ratio of FEV per second to total lung capacity established the presence of COPD) were shown as significant indicators of increased mortality. Specific occupations have increased the risk of lung diseases – Farmer’s Lung (caused by the increased inhalation of biologic dust particles) is one such disease. Braun et al. (Reference Braun, doPico, Tsiatis, Horvath, Dickie and Rankin1979) studies the long-term and clinical outcomes of this disease – this research (and other such research) coupled with our research could assist in competitive pricing of life insurance premiums and provide a framework for individual insurance pricing.

Covariance regression was originally proposed in Hoff & Niu (Reference Hoff and Niu2012). Algorithm 1 defines a modified version of the original algorithm where the local mean of the dependent variables can be estimated independently of the algorithm. This allows for increased flexibility in defining a local mean structure. A rank 1 covariance regression model seeks to calculate $ \boldsymbol{\Psi }$ , the base unattributable (or systemic or contemporaneous) covariance, and $ \mathbf{B}$ , the matrix of coefficients relating the covariance of $ \mathbf{Y}$ to the independent variables $ \mathbf{X}$ such that:

(30) \begin{equation} \text{cov}[\mathbf{Y}|\mathbf{X}] = \boldsymbol{\Psi } + \mathbf{B}\mathbf{X}\mathbf{X}^T\mathbf{B}^T. \end{equation}

In the following section, the independent mean-covariance regression model is introduced in translatable pseudocode in Algorithm 1. In this algorithm, Equation (31) is repeated from Hoff & Niu (Reference Hoff and Niu2012) which assumes the iterative inverses are well behaved. Within this software package, the pseudo-inverse is used to avoid unnecessary errors. The sections that follow this restatement of the algorithm (with a defined pseudocode algorithm and independent mean definition) are discussed with particular focus on variations of the Equation (31) for differing degrees of fitting and variable selection – these algorithmic variations are presented herein for the first time and applied to the same data as in Hoff & Niu (Reference Hoff and Niu2012) for direct comparison with the original technique both formally and graphically as in Figs. 10 and 11. These algorithmic variations are presented here for the reader to explore and develop – the implicit factor models developed in the previous section are ideal candidates for independent variables in the RCR framework owing to (most of) them yielding smoothed factors (compared to discontinuous explicit market factors such as the 10-year bond yield, etc.) for relatively smooth forecasted covariance transitions.

5.3 Independent mean-covariance regression

In Hoff & Niu (Reference Hoff and Niu2012), covariance regression is introduced along with the formal derivations. In this section, a modified and more robust version of the algorithm originally presented in Hoff & Niu (Reference Hoff and Niu2012) is presented and detailed in Algorithm 1.

In the absence of some stopping criterion, the modified covariance regression algorithm requires five inputs. Matrices $ \mathbf{A} \in \mathbb{R}^{(q\times{p})}$ and $ \mathbf{X}_1 \in \mathbb{R}^{({n}\times{q})}$ are the major modification to the original algorithm in that the means of the dependent variable, $ \mathbf{Y}$ , can be calculated and optimized independent of the variance calculation. This increased robustness is ideal for variance attributable to different frequency structures which are common in the financial setting. The mean matrix, $ \boldsymbol{\mu } \in \mathbb{R}^{(n\times{p})}$ is calculated as $ \boldsymbol{\mu } = \mathbf{X}_1{\mathbf{A}}$ with $q$ being the number of independent variables (implicit or explicit factors) used in the model, $p$ being the number of dependent variables, and $n$ being the number of observations.

With $ \mathbf{X}_2 \in \mathbb{R}^{(n\times{r})}$ being the matrix of independent variables which can be different to the basis matrix for mean construction, that is, different to $ \mathbf{X}_1$ , with $ \mathbf{Y} \in \mathbb{R}^{(n\times{p})}$ being the matrix of independent variables, $ \mathbf{B} \in \mathbb{R}^{({p}\times{n})}$ being the matrix of structural covariance coefficients, and with $ \mathbf{max\_iter}$ being the maximum allowed algorithmic iterations.

Algorithm 1 Covariance Regression

Using Algorithm 1, the data from Hoff & Niu (Reference Hoff and Niu2012), and the specification of the cubic B-spline knot points as in Hoff & Niu (Reference Hoff and Niu2012) (knots at $ 4, 11, 18$ ), one arrives at Equation (32) below with $ \mathbf{X}_2 = [\mathbf{1}, \textbf{age}^{\frac{1}{2}}, \textbf{age}]^T$ . Using these parameters, one can arrive at Figs. 8 and 9 to estimate the covariation of FEV and height as a function of age from Hoff & Niu (Reference Hoff and Niu2012):

(32) \begin{equation} \mathbf{B}_{direct} = \left [\begin{array}{rr} -2.60408& -10.45065\\ 1.43707& 6.44579\\ -0.14644& -0.82822 \end{array}\right ]. \end{equation}

Figure 8 Replication of Fig. 5 from (see p. 24, Hoff & Niu, Reference Hoff and Niu2012) plotting FEV data and variance in the left figure and height data and variance in the right figure as a function of age.

Figure 9 Replication of Fig. 6 from (see p. 25, Hoff & Niu, Reference Hoff and Niu2012) plotting FEV data and variance, height data and variance, and FEV and height data and correlation as a function of age.

5.4 LASSO covariance regression

Equation (31) is the direct calculation of $ \mathbf{B}$ , but in RCR $ \mathbf{B}$ can be calculated in several ways. $ \mathbf{B}_{LASSO}$ is calculated by minimizing the following:

(33) \begin{equation} \text{arg }\text{min}_{\mathbf{B}}\bigg ((\tilde{\mathbf{Y}} - \tilde{\mathbf{X}}_2\mathbf{B}^T)^T(\tilde{\mathbf{Y}} - \tilde{\mathbf{X}}_2\mathbf{B}^T) + \lambda _1{||\mathbf{B}||_1}\bigg ). \end{equation}

This form of optimization promotes sparsity as can be seen in Equation (34). In Equation (34), it can be noted that the variances and covariance of FEV and height concerning age are most dependent on age. As a result of the structure of Equation (30), the covariance matrix becomes a function of age squared:

(34) \begin{equation} \mathbf{B}_{\text{lasso}} = \left [\begin{array}{rr} -0.00000& -0.00000\\ -0.00000& -0.00000\\ -0.04537& -0.12452 \end{array}\right ]. \end{equation}

5.5 Ridge regression covariance regression

Ridge regression penalizes unjustifiably large coefficients in $ \mathbf{B}$ by optimizing the following:

(35) \begin{equation} \text{arg }\text{min}_{\mathbf{B}}\bigg ((\tilde{\mathbf{Y}} - \tilde{\mathbf{X}}_2\mathbf{B}^T)^T(\tilde{\mathbf{Y}} - \tilde{\mathbf{X}}_2\mathbf{B}^T) + \lambda _2{||\mathbf{B}||^2_2}\bigg ). \end{equation}

In this setting, by manipulating Equation (35), one can observe that ridge regression introduces a bias such that:

(36) \begin{equation} \text{arg }\text{min}_{\mathbf{B}}\bigg (||\tilde{\mathbf{Y}} - \tilde{\mathbf{X}}_2\mathbf{B}^T + \lambda _2\mathbf{B}||^2_2\bigg ). \end{equation}

By referring to Equation (31) and noting the structure of Equation (36), Equation (37) follows as a biased version of the original direct solution where $ \lambda _2$ can also be seen as a bias term rather than merely a penalty:

(37) \begin{equation} \mathbf{B} = \tilde{\mathbf{Y}}^T\tilde{\mathbf{X}}_2\Big (\tilde{\mathbf{X}}_2^T\tilde{\mathbf{X}}_2 + \lambda _2\mathbf{I}\Big )^{-1}. \end{equation}

By comparing Equation (32) against Equation (38), it can be noted that the coefficients have been greatly reduced as a result of the penalizing or biasing term:

(38) \begin{equation} \mathbf{B}_{\text{ridge}} = \left [\begin{array}{rr} 0.27839& 1.31838\\ 0.07562& -0.51146\\ -0.09374& -0.10844 \end{array}\right ]. \end{equation}

5.6 Elastic-net covariance regression

Elastic-net regression is a compromise between LASSO regression and Ridge regression where $ \mathbf{B}$ is optimized using the following:

(39) \begin{equation} \text{arg }\text{min}_{\mathbf{B}}\bigg ((\tilde{\mathbf{Y}} - \tilde{\mathbf{X}}_2\mathbf{B}^T)^T(\tilde{\mathbf{Y}} - \tilde{\mathbf{X}}_2\mathbf{B}^T) + \lambda _1{||\mathbf{B}||_1} + \lambda _2{||\mathbf{B}||^2_2}\bigg ). \end{equation}

This can be reparameterized (as is in Python) to:

(40) \begin{equation} \text{arg }\text{min}_{\mathbf{B}}\bigg ((\tilde{\mathbf{Y}} - \tilde{\mathbf{X}}_2\mathbf{B}^T)^T(\tilde{\mathbf{Y}} - \tilde{\mathbf{X}}_2\mathbf{B}^T) + \lambda _1r_1{||\mathbf{B}||_1} + \frac{1}{2}\lambda _1(1-r_1){||\mathbf{B}||^2_2}\bigg ), \end{equation}

with $ r_1$ being the ratio of the L1-norm to the L2-norm. Noting that Equation (39) can be viewed as a combination of Equations (33) and (36), it becomes clear why Equation (41) displays both sparsity and smaller parameters. This is seen as a favorable compromise between the two methods:

(41) \begin{equation} \mathbf{B}_{\text{elastic-net}} = \left [\begin{array}{rr} 0.00000& 0.00000\\ 0.02062& -0.04304\\ -0.05310& -0.12245 \end{array}\right ]. \end{equation}

5.7 Group-LASSO covariance regression

The group-LASSO regression promotes a sparse number of groups of parameters and in each group, there will be a sparse set of parameters:

(42) \begin{equation} \text{arg }\text{min}_{\mathbf{B}}\bigg ((\tilde{\mathbf{Y}} - \tilde{\mathbf{X}}_2\mathbf{B}^T)^T(\tilde{\mathbf{Y}} - \tilde{\mathbf{X}}_2\mathbf{B}^T) + \lambda _1{||\mathbf{B}||_1} + \lambda _g{\sum _{g\in \mathcal{G}}\sqrt{d_g}||B_g||_2}\bigg ). \end{equation}

In Equation (43), $ \lambda _g$ is the penalty term for the grouped coefficient terms, $ \mathcal{G}$ is a partition of the coefficients into subsets, $ d_g$ is the cardinality of each of the subsets (larger sets get penalized more than smaller sets), and $ B_g$ is the set coefficients within subgroup $ g$ . With a grouping of {0, 1, 1}, that is with the constant being in a group separate from age $^{\frac{1}{2}}$ and age, this results in the following parameter estimates:

(43) \begin{equation} \mathbf{B}_{\text{group-lasso}} = \left [\begin{array}{rr} 0.00823& -0.00091\\ 0.00577& -0.03660\\ -0.04977& -0.12262 \end{array}\right ]. \end{equation}

5.8 Subgradient covariance regression

Equation (44) can be seen as being identical to Equation (33) in Section 5.4. The differences lie in the method of approaching this solution. Subgradient descent or subgradient optimization allows optimization when functions are continuous, but non-differentiable:

(44) \begin{equation} \begin{aligned} \text{minimize } &||\mathbf{B}||_1\\ \text{subject to } &\tilde{\mathbf{X}}_2\mathbf{B}^T = \tilde{\mathbf{Y}}, \end{aligned} \end{equation}

To initialize the algorithm, several candidates can be used. In this algorithm, $ \mathbf{B}_0$ is initialized with Equation (31) which is the direct solution in Euclidean Norm Space. The next iteration of the solution, $ \mathbf{B}_{k+1}$ , is calculated using the following equation:

(45) \begin{equation} \mathbf{B}_{k+1} = \mathbf{B}_k - \alpha _k\bigg (\mathbf{I}_p - \tilde{\mathbf{X}}_2^T\big (\tilde{\mathbf{X}}_2\tilde{\mathbf{X}}_2^T\big )^{-1}\tilde{\mathbf{X}}_2\text{ sign}\big (\mathbf{B}_k\big )\bigg ), \end{equation}

with $ \text{sign}(\cdot )$ being defined as:

(46) \begin{equation} \text{sign}(x) = \left \{\begin{array}{rr} +1, & \text{ for } x \gt 0\\ 0, & \text{ for } x = 0\\ -1, & \text{ for } x \lt 0\\ \end{array}\right \}, \end{equation}

and with the incremental step size being a member of the square summable, but not summable family of step sizes. In this application, $ \alpha _k = \frac{10^{-2}}{k}$ with the formal conditions of the square summable step sizes required to satisfy the following conditions:

(47) \begin{equation} \alpha _k \gt 0\text{ }\forall \text{ }k \gt 0, \sum _{i=0}^\infty \alpha _k^2 \lt \infty, \text{ and } \sum _{i=0}^\infty \alpha _k = \infty . \end{equation}

With this technique applied to the dataset from Hoff & Niu (Reference Hoff and Niu2012) using $ 10$ iterations one arrives at:

(48) \begin{equation} \mathbf{B}_{\text{subgradient}} = \left [\begin{array}{rr} -1.77455& -4.37165\\ 1.44572& 3.53236\\ -0.28942& -0.71025 \end{array}\right ]. \end{equation}

To estimate the base unattributable (or systemic or contemporaneous) covariance, $ \boldsymbol{\Psi }$ , and the coefficients, $ \mathbf{B}$ , that relate the attributable (or structural) covariance to the independent variables one can apply the following code:

B_est, Psi_est = \

cov_reg_given_mean(A_est, basis, x, y, iterations, technique,

alpha, l1_ratio_or_reg, group_reg,

max_iter, groups)

Each of these techniques was applied to the original dataset from Hoff & Niu (Reference Hoff and Niu2012) and the results are plotted in Figs. 10 and11. In the code above, A_est and basis refer to the mean coefficients and basis to be used for means of the dependent variable, y, with x being the independent variable. The number of iterations of the core covariance regression algorithm is controlled by the iterations parameter. The technique is being used in estimating the $ \mathbf{B}$ matrix in each iteration of the total iterations. The options available for technique are “direct,” “lasso,” “ridge,” “elastic-net,” “group-lasso,” and “sub-gradient" which correspond to Sections 5.3, 5.4, 5.5, 5.6, 5.7, and 5.8, respectively. The unambiguous penalty term for LASSO, ridge, elastic-net, group-LASSO regression as well as the first value in the subgradient optimization is alpha.

The l1_ratio_or_reg controls the ratio of the l1-ratio to l2-ratio as in Equation (40) when performing elastic-net regression. The group-LASSO regression group penalty term is controlled by group_reg. While possibly misleading when compared against iterations, the max_iter parameter controls the number of iterations to be used in the subgradient descent algorithm as in Equation (45). The groups variable controls the grouping of the independent variables and the corresponding parameters using a vector of indices such as in Section 5.7. Finally, the test_lasso is a Boolean variable that controls whether the alpha value for LASSO regression is optimized.

Figure 10 Replication of Fig. 5 from (see p. 24, Hoff & Niu, Reference Hoff and Niu2012) plotting FEV data and variance in the left figure and height data and variance in the right figure as a function of age with RCR alternate estimate.

Figure 11 Replication of Fig. 6 from (see p. 25, Hoff & Niu, Reference Hoff and Niu2012) plotting FEV data and variance, height data and variance, and FEV and height data and correlation as a function of age with RCR alternate estimate.

5.9 Closing remarks on lung capacity versus physiology as a cost-saving benefit to insurers

Section 5.1 notes the decrease in costs to insurers from the early detection of a number of degenerative lung diseases. Several of the more well-known conditions that affect human lungs are listed in the introduction to Section 5. The physiological variable that is perceived as the most highly correlated with most of these conditions is a person’s age. This research has provided a framework for the early detection of the presence of a number of degenerative lung diseases using FEV measurement (as outlined in Section 5.2) as well as other relevant variables, besides age. With this framework, one can use any number of lifestyle or physiological variables to assess the presence of outliers in FEV measurements which can lead to early detection and significantly reduce the costs to insurers through comparatively cheaper early-stage treatment compared against late-stage treatment or end-of-life therapy.

6.1 Unfunded and funded pension schemes

Chapman et al. (Reference Chapman, Gordon and Speed2001) notes the lower and more certain returns of bond investment (as a proxy for fully unfunded) versus the higher and less certain returns of equity investment (as a fully funded pension proxy). Chapman et al. (Reference Chapman, Gordon and Speed2001) goes further in noting that employees, the enterprise, the stakeholders (shareholders), and the government seldom benefit from the same schemes. The stakeholders in the company (and the associated share price) should be taken into account when setting the contribution rates for the unfunded (or defined benefits) portion of the pension as well as the investment strategies of the companies once funds have been attained. Chapman et al. (Reference Chapman, Gordon and Speed2001) also notes the need for clarity for the employees and company as the funded portion of the pensions are secured versus the relatively unsecured portion of the pension that is unfunded should there be unforeseen market shocks such as the 2008–2009 financial crisis or the SARS-CoV-2 pandemic.

Splinter (Reference Splinter2017) notes the broader economic benefits of maintaining an at least partially unfunded pension scheme for state employees in which it provides a stabilizing effect on a nation’s economy by providing governments with funds without needing to either raise taxes or cut benefits. There is a downside risk in that the government (and by extensions the employees) then relies on the same assumptions of unfunded pensions in that the economy will continue to grow and the money can be provided or returned at a later date. Splinter (Reference Splinter2017) notes that the main contributor to the under-funding of the public employees pensions is the unfunded portions as a result of an increase in state liabilities with an insufficient associated contribution by employees and the government. Whereas, the funded portions properly invested in diversified equity portfolios’ returns were in excess of expectations. Other causes were noted, but the large accumulation of liabilities was seen as the main contributor to less-than-favorable defined benefit pension funding.

Grubbs (Reference Grubbs1999) notes the difficulty experienced by some private firms in adequately funding defined benefits schemes before adequate legislation and planning was enforced. These failed fully defined benefit plans (or unfunded pensions) relies on inadequate contributions from employees and the companies. These resulted in a number of disastrous instances such as companies canceling all their unfunded or defined benefits schemes or only being able to continue paying a small portion of the retirees at the expense of the future retirees. An extended study on this topic was conducted in Trowbridge (Reference Trowbridge1952).

With the problems with pure funding and unfunded pensions noted above and with the rise of demographic mismatches in developed populations discussed in the next section, Section 6.2, Devolder & de Valeriola (Reference Devolder and de Valeriola2019) notes that neither funded (referred to as defined contributions or DC) nor unfunded (referred to as defined benefits or DB) are adequate as risk-sharing compromises between generations. This, among other things, motivated the rise of hybrid strategies. The optimization in Devolder & de Valeriola (Reference Devolder and de Valeriola2019) of the ratio of funded to unfunded schemes is considered from the perspective of minimizing the disparity between retirement living conditions of different generations rather than some risk-returns trade-off of individual investors which is performed in this case study.

Wang & Lu (Reference Wang and Lu2019) studies hybrid pensions strategies to address the problem of risk-sharing between generations from a stochastic perspective with adjusted contributions and benefits depending on the performances of the pensions. Using the stochastic optimal control approach, closed-form solutions are derived using both a quadratic loss function and an exponential loss function. Various market parameters are adjusted with the contributions and benefits adjusted accordingly to demonstrate the effectiveness of risk-sharing between different generations when compared against purely funded or unfunded plans. Some simplifying assumptions were made in the analysis (such as assuming mortality rates remained constant over time) to allow the calculation of closed-form solutions.

Further, Blommestein et al. (Reference Blommestein, Janssen, Kortleve and Yermo2009) notes that hybrid pension strategies perform better for individual private investors in terms of funding ratios and replacement ratios. The simulations performed by Blommestein et al. (Reference Blommestein, Janssen, Kortleve and Yermo2009) show that hybrid strategies perform best in terms of sustainable financial risk-sharing between generations. It is found that in relatively stable environments where over-funding of hybrid strategies takes place, the most effective way of sharing risk is through conditional index-linked investing. It is shown how hybrid plans offer an agreeable compromise between DB plans with no benefit risk to the individual, but risk to the employer, and DC plans where the individual bears the entirety of the risk through investment decisions, inflation, and increased longevity.

Hoevenaars & Ponds (Reference Hoevenaars and Ponds2008) evaluates these hybrid options from the balance sheet perspective as the sum of embedded generational options. This leads to the conclusion that any policy change of an individual inevitably leads to the transfer of risk between generations. It is found, maybe not surprisingly, that shifts in the investment pools of pensions to less risky assets shifts the risks from older generations to younger generations. Whereas, a shift in weightings between funded and unfunded within the hybrid strategy from flexible contributions and fixed benefits to fixed contributions and variable benefits shifts the risk from younger generations to older generations.

6.2 Hybrid funded–unfunded UK pension schemes to address demographic changes

The changing demographics of the more-developed countries in the world (proportionally larger aging working populations than developing countries) are putting significant strain on unfunded or pay-as-you-go (PAYG) pension systems – see Miles (Reference Miles2000). The motivation and evidence for moving away from PAYG pensions is based on developed economies having efficient market portfolios that significantly exceed GDP growth and funded pensions lead to returns that are more robust to temporal distortions in the availability of working-age people within the population. These reasons, while arguably logical, are not in themselves conclusive for transfers from unfunded to funded pension systems for more developed nations.

Miles (Reference Miles1998) notes that for this shift from unfunded to funded pensions, the population would need more savings. It is also noted that given the aging populations and the unfunded pensions on which they depend, significant leveraged funding would be required by the governments to fund these shifts in policy or there would naturally be some losers when shifting systems owing to the funding needing to come from somewhere. It is, however, noted in Miles (Reference Miles1998) that the UK is further along in this process than most of its mainland Europe counterparts which face significant difficulties that have no doubt been exacerbated by the 2007–2008 financial crisis and SARS-CoV-2 pandemic-induced recessions since the its publication. These demographic changes and the resulting difficulties they cause for governments in more developed countries to shift from unfunded to funded pension schemes have given rise to significant literature on the topic such as Disney (Reference Disney1996), Feldstein (Reference Feldstein1996), Kotlikoff (Reference Kotlikoff1996), Mitchell & Zeldes (Reference Mitchell and Zeldes1996), Roseveare et al. (Reference Roseveare, Leibfritz, Fore and Wurzel1996), Feldstein & Samwick (Reference Feldstein, Samwick and Feldstein1998), and Miles et al. (Reference Miles, Timmermann, de Haan and Pagano1999).

Sinn (Reference Sinn1999) supports the full or partial shift away from PAYG pension systems as a way or addressing the growing demographic-related problems in more developed countries by replacing the loss of labor with capital. These demographic issues (larger ratio of retirees to working-age population, the associated costs of this relative aging population, and the drastically falling birth rates) put strain on PAYG systems in a similar manner as pyramids being dependent on a broader base.

Martell et al. (Reference Martell, Kioko and Moldogaziev2013) also notes the difficult decision facing governments when funding pensions (by matching personal contributions) using debt increases the scheme’s funding ratio which leads to a country’s credit rating and global perception decreasing. The majority of the dangers of under-funding pensions for governments are related to the costs associated with an aging population. The costs are increasing as a result of the increasing ratio of retirees to the working population, the increasing life expectancy of retirees, and the increasing healthcare costs associated with an aging population living longer.

6.3 Equal RPP weighting strategy

RPP portfolio weighting strategies, proposed in Maillard et al. (Reference Maillard, Roncalli and Teïletche2010), have had empirical successes when back-testing on equities over previous financial crises. The RPP weights are calculated as follows:

(49) \begin{equation} \boldsymbol{\omega }_{\text{med}} = \text{argmin}\sum _{j=1}^N\Bigg ({\omega _j}\frac{(\mathbf{\Sigma }_{\text{med}}\boldsymbol{\omega })_j}{\boldsymbol{\omega }^T\boldsymbol{\Sigma }_{\text{med}}\boldsymbol{\omega }} - b_j\Bigg ), \end{equation}

with $ \boldsymbol{\omega }_{\text{med}}$ being the RPP weighting vector, $ \mathbf{\Sigma }_{\text{med}}$ being the median of the forecasted covariance, $ N$ being the number of assets, and $ b_j$ being the relative proportion of risk allocated to asset $j$ . Only the equal RPP portfolio weighting strategy ( $ b_j = \frac{1}{N}$ ) is explored in this paper and package, with the other constraints being:

(50) \begin{equation} \sum _{j=1}^{N}\omega _j=1, \text{ and} \end{equation}

(51) \begin{equation} \sum _{j=1}^{N}\omega _j\boldsymbol{I}_{\{\omega _j\geq{0}\}}\geq{k_{\text{long}}} \text{ } \text{ for } k_{\text{long}} \geq 1.3, \text{ and} \end{equation}

(52) \begin{equation} \sum _{j=1}^{N}\omega _j\boldsymbol{I}_{\{\omega _j\leq{0}\}}\geq{k_{\text{short}}} \text{ } \text{ for } k_{\text{short}} \leq -0.3. \end{equation}

The final two constraints are referred to as long/short equity constraints that allow a certain level of shorting ( $ k_{\text{short}}$ ) to finance further long positions. These two weighting strategies are not exhaustive and are included here for context, clarity of case studies in the package, and completeness.

6.4 Mean-variance utility

We follow the model in Dutta et al. (Reference Dutta, Kapur and Orszag2000) with the stochastic return of the portfolio over time increment $ t - 1$ to $ t$ being $ P_t = w_tR_{f,t} + (1 - w_t)R_{u,t}$ with $ w_t$ being the proportion of the portfolio in the diversified RPP portfolio based on either high-, mid-, or low-frequency structures with the stochastic returns of said portfolio being $ R_{f,t}$ and $ R_{u,t}$ being the stochastic returns of the unfunded portion which we approximate with the quarterly UK GDP. The resulting mean-variance utility function is such that:

(53) \begin{equation} \mathbb{E}[U(P_t)] = \mathbb{E}[P_t] - \frac{\gamma }{2}\mathbb{V}ar(P_t). \end{equation}

with $ \gamma$ being the risk aversion level to be discussed later. The resulting expected returns and variance of the portfolio are as follows:

(54) \begin{equation} \begin{aligned} \mathbb{E}[P_t] &= w_t\mathbb{E}[R_{f,t}] + (1 - w_t)\mathbb{E}[R_{u,t}]\\ &= w_t\mu _{f,t} + (1 - w_t)\mu _{u,t}, \end{aligned} \end{equation}

and

(55) \begin{equation} \mathbb{V}ar(P_t) = w_t^2\sigma ^2_{f,t} + (1 - w_t)^2\sigma ^2_{u,t} + w_t(1 - w_t)\sigma _{fu,t}, \end{equation}

respectively, with the realized returns of the funded portfolio over time increment $ t - 1$ to $ t$ being $ \mu _{f,t}$ (for either the high-, mid, or low-frequency covariance forecasts), the realized returns of the unfunded portfolio over the same increment being $ \mu _{u,t}$ , $ \sigma ^2_{f,t}$ being the variance of the funded portfolio, $ \sigma ^2_{u,t}$ being the variance of the unfunded portion, and $ \sigma _{fu,t}$ being the covariance between the two returns. By taking the derivative of Equation (53) with respect to $ w_t$ and setting to zero, it can be shown that:

(56) \begin{equation} w_t = \frac{\mu _{f,t} - \mu _{u,t} + \gamma (\sigma ^2_{u,t} - \sigma _{fu,t})}{\gamma (\sigma ^2_{f,t} + \sigma ^2_{u,t} - 2\sigma _{fu,t})}. \end{equation}

6.5 Different risk aversions in FTSE 100 hybrid pension schemes

In this study, we take the daily logarithmic returns of the constituents of the FTSE 100 from 1 October 2003 until 1 April 2023 (where available) to construct RPP portfolios using the high-, mid-, and low-frequency structures to forecast the covariance of the constituents and weight the portfolios based on these forecasts. The high-, mid-, and low-frequency structures from the daily logarithmic returns isolated from 1 October 2003 until 31 December 2003 are regressed against the daily logarithmic returns for the same time series from 1 January 2004 until 31 March 2004.

The high-, mid-, and low-frequency structures isolated from the daily logarithmic returns isolated from 1 January 2004 until 31 March 2004 are then used to forecast the covariance of the constituents from 1 April 2004 until 30 June 2004. The portfolio is then weighted based on the median of these daily covariance forecasts using the different frequency structures. The case study and resulting returns therefore begin from 1 April 2004 in Fig. 12. This process is repeated until the available GDP proxy of the FTSE 100 quarterly returns is exhausted which is 1 April 2023.

Rather surprisingly (considering the RPP optimization focusing on the minimization of risk), the RPP portfolios weighted based on the high-, mid-, and low-frequency forecasts of the upcoming investment period’s covariance structure have higher returns when compared against the FTSE 100 and less desirable variance and value-at-risk (VaR). The cumulative returns of the portfolios, besides the FTSE 100 proxy, are much more susceptible to market-wide shocks such as those seen during the 2007–2008 financial crisis and the SARS-CoV-2 pandemic-induced recession. Based on the variance, the high-frequency content is more accurate when used to forecast the covariance than the mid-frequency and low-frequency content, but the VaR indicates that the low-frequency content is more susceptible to extreme shocks.

Figure 12 Plot of the cumulative returns of funded and unfunded pensions from 1 April 2004 until 1 April 2023.

Noting the risk aversion level ( $\gamma$ ) in Equations (53) and (56), we see that as the level of risk aversion is increased sequentially from $1.2$ in Table 2 to $1.4$ in Table 3 and, finally, to $1.6$ in Table 4 the weightings in the respective funded pension portions of the hybrid schemes decrease. This study indicates that a desirable compromise between the level of acceptable risk and the expected returns of the pensions can be attained by weighting the funded portions of the pensions using RPP based on covariance forecasts and increasing or decreasing the risk aversion level to meet an individual’s need regarding the proportion of the pension scheme to be funded or unfunded.

Table 2. Summary statistics for quarterly returns of hybrid pensions schemes using risk aversion level of $\gamma = 1.2$

Table 3. Summary statistics for quarterly returns of hybrid pensions schemes using risk aversion level of $\gamma = 1.4$

Table 4. Summary statistics for quarterly returns of hybrid pensions schemes using risk aversion level of $\gamma = 1.6$

6.6 Closing remarks on hybrid pensions schemes

In Section 6.1, we note works that are critical of both purely funded and purely unfunded pensions schemes. In general (depending on the prevailing economic conditions over the lifetime of the scheme), the funded schemes have higher and less certain returns with the returns of the unfunded schemes being lower and more certain. Numerous works cited in this section allude to the necessity for some desirable compromise between these two schemes. Section 6.2 notes not only that hybrid pension schemes would benefit the holders of these schemes, but they would also address arising demographic issues in the majority of more-developed countries as the shrinking workforces and growing retired populations are putting siginificant financial strain on traditionally unfunded pension schemes.

We introduce the RPP weighting strategy in Section 6.3 which we use in this case study to allocate portfolio weights. The mean-variance utility methodology of hybrid pension schemes which we use to allocate the ratio of our unfunded portion (using the FTSE 100 as a proxy) and our funded portion (using difference frequency covariance forecasts) is outlined in Section 6.4. The construction of the unfunded proxy and the funded portfolios using different frequency forecasts is given exposition in Section 6.5. We display the results of both purely unfunded and purely funded pension schemes using mean returns as well as various common portfolio risk measures in Table 1.

We note that the funded portfolios have higher returns as well as significantly higher variances, VaRs, and CVaRs. This leads to the hybrid pension scheme results in the following tables. By noting Tables 2, 3, and 4, one can infer that the high-frequency structures used to forecast the forthcoming investment period’s covariance structure are more responsive to changes in the covariance structures and therefore require a lower relative weighting (when compared against using mid- and low-frequency structures) in the funded portions of the portfolios to lower the risk measures such portfolio variance, VaR, and CVar. The mid-frequency structures are less responsive than the high-frequency structures while being more responsive than the low-frequency and constant structures. This forms a hierarchy of performance in both overall returns and the risk measures.