An Empirical Implementation of the Shadow Riskless Rate

Abstract

We address the problem of asset pricing in a market where there are no risky assets. Previous work developed a theoretical model for a shadow riskless rate (SRR) for such a market, based on the drift component of the state-price deflator for that asset universe. Assuming that asset prices are modeled by correlated geometric Brownian motion, in this work, we develop a computational approach to estimate the SRR from empirical datasets. The approach employs principal component analysis to model the effects of individual Brownian motions, singular value decomposition to capture abrupt changes in the condition number of the linear system whose solution provides the SRR values, and regularization to control the rate of change of the condition number. Among other uses such as option pricing and developing a term structure of interest rates, the SRR can be used as an investment discriminator between different asset classes. We apply this computational procedure to markets consisting of various groups of stocks, encompassing different asset types and numbers. The theoretical and computational analysis provides the drift as well as the total volatility of the state-price deflator. We investigate the time trajectory of these two descriptive components of the state-price deflator for the empirical datasets.

1. Introduction

One definition of a riskless (safe) asset is “an asset that is (almost always) valued at face value without expensive and prolonged analysis” (Gorton 2017). A safe asset changes hands with “no questions asked” (Holmstrom 2015). Riskless assets play a crucial role in financial markets, both in practice and in financial theory. Safe assets are used by financial entities to satisfy regulatory requirements, serve as pricing benchmarks, act as collateral in financial transactions, and aid in the development of asset and derivative pricing. Gorton et al. (2012) have shown that, while the total volume of US financial assets has increased 250% since 1952, the percentage of safe debt in the US economy has remained relatively stable at approximately 32%. This stability in the percentage of safe debt obscures the fact that the composition of safe assets has shifted from consisting primarily of government debt (US Treasury securities) and cash (demand deposits) to an ever-increasing reliance on innovative financial instruments (e.g., money market mutual fund shares, commercial paper, repurchase agreements, and securitized debt such as asset-backed and mortgage-backed securities) developed by the shadow banking system. The growth in the volume of financial assets, outpacing the availability of traditional safe assets, reflects the fact that shadow institutions are not subject to the same regulations as depository banks, do not need to keep the same level of financial reserves relative to market exposure, and can maintain higher levels of financial leverage.

This vulnerability became apparent during the 2008 global financial crisis. A consequence of subsequent regulations, such as the US Dodd–Frank Act and the BCBS Basel III Accord, has been a reduction in the availability of safe assets. However, a number of studies (Aoki et al. 2014; Bocola and Lorenzoni 2023; Caballero and Simsek 2020, 2021; Caballero and Farhi 2013; Caramp 2024; Eggertsson and Krugman 2012; Gorton et al. 2022; Gourinchas and Jeanne 2012) have focused on the challenges that a shortage of safe assets could pose for the stability of the financial system and monetary policy. The findings of these researchers suggest that the imposition of a liquidity coverage ratio, which requires banks to back their short-term debt dollar-for-dollar with treasury securities, risks both lowering interest rates and making financial panics more frequent. Caballero and Farhi (2018) described how a self-reinforcing demand for safe assets that are in short supply could prolong a recession into outright stagnation. Altinoglu (2023) introduced a theory of safe asset creation to study the interaction between systemic risk, including financial leveraging, and aggregate demand, particularly how the restriction of the safe asset supply through external constraints affects these dynamics.

As early as 1972, economists considered equilibrium markets to have no riskless assets. For example, Black (1972) showed that the capital asset pricing model is still valid even without a riskless asset. Recognizing that the continued development of non-traditional safe assets is a “train in motion that probably can’t be stopped”, it is necessary to replace qualitative definitions of a safe asset with a strong theoretical (and consequent empirical) construct. In response to this, Rachev et al. (2017) developed a riskless rate and a corresponding riskless asset for a market consisting solely of N risky assets. Theoretically, this riskless rate is (the negative of) the drift term of the state-price deflator for that market. The authors first considered the case where the price dynamics of the N assets were correlated geometric Brownian motions (GBMs) driven by $N-1$ Brownian motions. In addition to this continuous diffusion case, they also explored scenarios in which the price dynamics of the assets were determined by correlated jump-diffusion processes, diffusion with stochastic volatilities, and geometric fractional Brownian and Rosenblatt processes. For each chosen model of price dynamics, they derived the Black–Scholes–Merton equations to price a European contingent claim, constructing a tradable, perpetual derivative that serves as a proxy for a riskless asset whose price evolves according to this riskless rate. The resultant market, consisting of the N risky assets and this perpetual derivative, is complete and arbitrage-free. As the riskless rate derived for this perpetual derivative is applicable to the selected market class, and as it reflects the efforts of shadow bank portfolio managers to create riskless rates using their available assets, it is referred to as the shadow riskless rate (SRR).

To date, there has been no empirical investigation associated with this theoretical development. This paper aims to develop a practical numerical implementation of the approach proposed by Rachev et al. (2017) under the correlated GBM assumption (the simplest model under which to develop, test, and improve the needed numerical methods). For N assets whose prices are driven by $N-1$ Brownian motions, the theoretical method for deriving an SRR, as introduced by Rachev et al. (2017), is summarized in Section 2.1. Theoretically, the SRR is obtained as the unique solution to a linear system. In Section 2.2, given a window of empirical historical returns for the assets, we describe how to use principal component analysis to model the actions of the Brownian motions. Based on the PCA analysis, we derive two methods for generating the components of the coefficient matrix required for that linear system. In Section 2.3, we develop the moving window method to generate a time series of SRR values. We demonstrate numerically unstable behavior in the time series when the coefficient matrix exhibits rapid changes in the condition number. In Section 2.4 and Section 2.5, we develop an approach to provide a regularized solution when significant changes in the condition number of the linear system are detected. Section 2.6 is devoted to comparing solutions obtained using two different methods for generating the coefficient matrix. As the SRR is based on the drift component of a state-price deflator, an additional outcome of this analysis is the estimation of the total volatility of this deflator. The numerical estimation of this is discussed in Section 2.7.

The application of the numerical method to empirical datasets is described in Section 3. Section 3.1 computes estimates of the SRR for four different groups of N stocks as well as a group of N exchange-traded funds. Section 3.2 addresses the effect of the number, N, of stocks in the group. Section 3.3 focuses on the behavior of the state-price deflator, plotting the trajectory over time of its drift and volatility. Section 4 concludes the paper with a discussion of future research directions and a proposal for an alternate empirical definition of a riskless rate and its proxy asset.

2. Method

2.1. Derivation of the Shadow Riskless Rate

Consider a portfolio consisting of N risky assets ${\mathcal{S}}^{\left(j\right)}$, where $j=1,\dots ,N$, $N\ge 2$. Asset ${\mathcal{S}}^{\left(j\right)}$ has a price $S_{jt}$, $t\ge 0$, with its dynamics driven by $N-1$ Brownian motions $W_{kt}$, where $k=1,\dots ,N-1$; the dynamics are governed by the following stochastic differential equation:

(1)$d{S}_{jt}={\mu }_j{S}_{jt}dt+S_{jt}\sum _{k=1}^{N-1}{\sigma }_{jk}d{W}_{kt},t\ge 0,{\mu }_j>0,{\sigma }_j>0,j=1,\dots ,N,$

with the initial conditions $S_{j{t}_0}=S_{j0}>0$, $j=1,\dots ,N$. We define the mean return vector $\mu =[{\mu }_1,\dots ,{\mu }_N]$ and the variance–covariance matrix $\Sigma {\Sigma }^T$, where ${\Sigma }_{jk}={\sigma }_{jk}$, $j=1,\dots ,N$, $k=1,\dots ,N-1.$ The Brownian motions generate a filtered probability space $(\Omega ,\mathcal{F}$, $\mathbb{F}=\{{\mathcal{F}}_t,t\ge 0\},\mathbb{P})$ representing the natural world. The market of risky assets ${\mathcal{S}}^{\left(j\right)}$, $j=1,\dots ,N$, with price dynamics described by (1), will be complete if and only if there exists a unique state-price deflator ${\pi }_t$, $t\ge 0$, on $\mathbb{P}$ (Duffie 2001, sct. 6D) with dynamics given by the Itô process, as follows:

(2)$d{\pi }_t={\mu }_{\pi }{\pi }_{t}dt+{\pi }_t\sum _{k=1}^{N-1}{\sigma }_{\pi k}d{W}_{kt},t\ge 0.$

The existence and uniqueness of ${\pi }_t$ is equivalent to the statement that each deflated process $S_{jt}{\pi }_t$, $j=1,\dots ,N$, is a $\mathbb{P}-$martingale. This leads to the requirement that the linear system:

(3)${\mu }_j+{\mu }_{\pi }+\sum _{k=1}^{N-1}{\sigma }_{jk}{\sigma }_{\pi k}=0,j=1,\dots ,N,$

have unique solutions for ${\mu }_{\pi }$ and ${\sigma }_{\pi k}$, $k=1,\dots ,N-1$. In matrix notation (3) can be written as follows:

(4)$\Phi x=\mu ,$

where $\mu$ is the column vector ${[{\mu }_1,\dots ,{\mu }_N]}^{\mathrm{T}}$, x is the solution vector ${[-{\mu }_{\pi },{\sigma }_{\pi 1},\dots ,{\sigma }_{\pi (N-1)}]}^{\mathrm{T}}$, and the matrix $\Phi =[11_N,-\Sigma ]$, with $11_N$ being the $(N\times 1)$ column vector composed of ones. The SRR $\nu$ is given by $\nu =-{\mu }_{\pi }$ (Rachev et al. 2017). From (4), it has the analytic solution

(5)$\nu ={\displaystyle \frac{det{\Phi }_{\mu }}{det\Phi }},$

where ${\Phi }_{\mu }=[\mu ,-\Sigma ]$. Rachev et al. (2017) constructed the price process $B_t$ of a tradable, perpetual European contingency contract $\mathcal{B}$ that serves as a riskless bond in the arbitrage-free, complete market $({\mathcal{S}}_1,\dots ,{\mathcal{S}}_N,\mathcal{B})$. The price dynamics of $\mathcal{B}$ obeys the following:

(6)$d{B}_t=\nu B_{t}dt.$

The analytic solution for each standard deviation ${\sigma }_{\pi k}$, $k=1,\dots ,N-1$, can be expressed similarly as a ratio of determinants. Thus, while ${\mu }_{\pi }=-\nu$ measures the drift component of the state-price deflator, the vector $[{\sigma }_{\pi 1},\dots ,{\sigma }_{\pi (N-1)}]$ estimates the total variance ${\sigma }_{\pi }^2$ of the state-price deflator, as follows:

(7)${\sigma }_{\pi }^2=\sum _{k=1}^{N-1}{\left({\sigma }_{\pi k}\right)}^2.$

We will refer to ${\sigma }_{\pi }=\sqrt{{\sigma }_{\pi }^2}$ as the (total) volatility of the state-price deflator.

We will make some critical observations of the SRR and the total volatility of the state-price deflator by briefly describing the solution when $N=2$. (See Rachev et al. (2017)). Equation (1) are as follows:

(8)$\begin{array}{cc}\hfill d{S}_{1t}& ={\mu }_1{S}_{1t}dt+{\sigma }_1{S}_{1t}d{W}_t,{\mu }_1>0,{\sigma }_1>0,t\ge 0,\hfill \\ \hfill d{S}_{2t}& ={\mu }_2{S}_{2t}dt+{\sigma }_2{S}_{2t}d{W}_t,{\mu }_2>0,{\sigma }_2>0,t\ge 0;\hfill \end{array}$

Equation (2) is as follows:

(9)$d{\pi }_t={\mu }_{\pi }{\pi }_{t}dt+{\sigma }_{\pi }{\pi }_{t}d{W}_t,t\ge 0;$

and the linear system (3) can be written as follows:

(10)$-{\mu }_{\pi }={\mu }_1+{\sigma }_1{\sigma }_{\pi }={\mu }_2+{\sigma }_2{\sigma }_{\pi },$

leading to the following solutions:

(11a)$\begin{array}{cc}\hfill \nu =-{\mu }_{\pi }& ={\displaystyle \frac{{\mu }_1{\sigma }_2-{\mu }_2{\sigma }_1}{{\sigma }_2-{\sigma }_1}},\hfill \end{array}$

(11b)$\begin{array}{cc}\hfill {\sigma }_{\pi }& ={\displaystyle \frac{{\mu }_1-{\mu }_2}{{\sigma }_2-{\sigma }_1}}.\hfill \end{array}$

Equation (10) provides a useful identity between the market prices of risk for two assets, as follows:

(12)${\displaystyle \frac{{\mu }_1-\nu }{{\sigma }_1}}=-{\sigma }_{\pi }={\displaystyle \frac{{\mu }_2-\nu }{{\sigma }_2}},$

reinforcing the view of $\nu$ as a riskless rate. Without loss of generality, we can assume $\Delta \sigma ={\sigma }_2-{\sigma }_1>0$. Defining $\Delta \mu ={\mu }_2-{\mu }_1$, (11b) and (11a) can be written as follows:

(13)${\sigma }_{\pi }=-{\displaystyle \frac{\Delta \mu }{\Delta \sigma }},\nu =-{\mu }_{\pi }={\mu }_1\left(1+{\displaystyle \frac{{\sigma }_{\pi }}{{\mu }_1/{\sigma }_1}}\right).$

As ${\mu }_1>0$, ${\sigma }_1>0$, and $\Delta \sigma >0$, we see that ${\sigma }_{\pi }$ can either be positive or negative depending on the sign of $\Delta \mu$. If ${\sigma }_{\pi }$ is positive, then $\nu$ is positive (and ${\mu }_{\pi }$ is negative). If ${\sigma }_{\pi }$ is negative, $\nu$ can be either positive or negative depending on the magnitude of the ratio ${\sigma }_{\pi }/({\mu }_1/{\sigma }_1)$. For $N>2$ risky assets, the determination of the sign of the solution (5) is more complicated but can be either positive or negative.

2.2. Calibration to Historical Data

In practice, the estimation of the SRR via (4) involves the estimation of the components ${\mu }_j$ and ${\sigma }_{jk}$ of $\mu$ and $\Sigma$ in the discrete setting $dt\approx \Delta t$. In discrete time, the log return $r_j$ for asset j will have an instantaneous mean ${\mu }_j$ and instantaneous standard deviation vector $[{\sigma }_{j1},\dots ,{\sigma }_{j(N-1)}]$. Estimating these components requires a method to retrieve the $N-1$ Brownian motions that drive the uncertainty in the return dynamics. As the distributional assumption embodied in (1) is multivariate geometric Gaussian, we use principal component analysis (PCA) to do this. A brief summary of PCA is provided in Appendix A. Here, we employ the matrix row and column vector notation introduced in Appendix A.

Let R denote the $M\times N$ matrix of daily log returns $r_{mj}$ realized by the N risky assets over M historical trading days (i.e., $\Delta t=1$). We estimate the mean log-return value ${\mu }_j$ of each asset j from the historical data. The historical log returns can be adjusted by subtracting their historical mean values, resulting in a modified matrix R, where each column has a zero mean.

Applying PCA to R results in N-ordered eigenvector–eigenvalue pairs. As the multivariate price process is driven by the multivariate geometric Brownian motion, the returns $r_{mj}$ are distributed in a hyper-ellipsoid in ${\mathbb{R}}^N$. The column eigenvectors $w_{\circ 1},\dots ,w_{\circ N}$ from the PCA, arranged in order of the variance they explain, represent the principal axes of this hyper-ellipsoid. The components $p_{mj}$, $m=1,\dots ,M$ of the principal component column vector $p_{\circ j}=R{w}_{\circ j}$ describe the marginal distribution of the vector of correlated returns $(r_{m1},\dots ,r_{mN})$ on the j’th principal axis of the hyper-ellipsoid. As there are $N-1$ Brownian motions driving the return distribution, we make the approximation that the specific actions of these $N-1$ Brownian motions on the observed dataset are captured by the principal components along the first $N-1$ principal axes, $w_{\circ j}$, $j=1,\dots ,N-1$. (We ignore the distribution along $w_{\circ N}$, as the principal axes explain the smallest variance of the observed dataset). Based on (A2), we have the following approximation:

(14)${\sigma }_{jk}=\sqrt{{\lambda }_k}{w}_{jk},j=1,\dots ,N,k=1,\dots ,N-1.$

There is a second method to approximate the values ${\sigma }_{jk}$. Again, assuming that the actions of the $N-1$ Brownian motions on the observed dataset are described by the distributions of the principal components observed along the first $N-1$ principal axes, we can approximate the following:

(15)$r_{mj}-{\mathrm{E}}_m\left[r_{mj}\right]=\sum _{k=1}^{N-1}{\sigma }_{jk}{\overline{P}}_{ik},m=1,\dots ,M,j=1,\dots ,N,$

where:

(16)${\overline{P}}_{ik}={\displaystyle \frac{P_{ik}-{\mathrm{E}}_i\left[P_{ik}\right]}{\sqrt{{\mathrm{Var}}_i\left[P_{ik}\right]}}},$

with ${\mathrm{E}}_m\left[r_{mj}\right]$ and ${\mathrm{E}}_i\left[P_{ik}\right]$ denoting the average of the indicated column, and ${\mathrm{Var}}_i\left[P_{ik}\right]$ denoting the column variance. Solutions to the values ${\sigma }_{jk}$ in (15) are obtained via linear regression.¹

2.3. Computation of a Time Series for the SRR

The calibration described in Section 2.2 requires a historical window (of size M days) to compute a single value of the SRR $\nu$ via (4). Using a standard moving window procedure, we can develop a time series of historical SSR values. Given a group of N assets, we compute a value for ${\nu }_t$ for the date t via the following steps.

1. Assemble the log-return matrix $R_t$ for these N assets over the historical window $\{t-(M-1),\dots ,t-1,t\}$ of M trading days.

2. Estimate the vector ${\mu }_t$ from the historical data and subtract the respective mean values from each column of the matrix $R_t$.

3. Perform a PCA of $R_t$, producing the ordered eigenvalue and eigenvector pairs, ${\lambda }_j^{\left(t\right)}$, $w_{\circ j}^{\left(t\right)}$, $j=1,\dots ,N$.

4. Generate the matrix ${\Sigma }_t$ composed of column vectors ${\sigma }_{\circ k}^{\left(t\right)}$, $k=1,\dots ,N-1$ using either (4.i) the direct solution (14) or (4.ii) the regression solution to (15).

5. Form ${\Phi }_t$.

6. While ${\nu }_t$ can be obtained analytically from (5) and the standard deviations ${\sigma }_{\pi k}$ can be obtained from analogous determinant ratios, it is more numerically efficient and stable to solve (4) using LU decomposition with pivoting.

Steps 1 to 6 generate the time series, as follows:

(17)${\nu }_t=x_{1,t},{\sigma }_{\pi k,t}=x_{k,t},k=2,\dots ,N.$

Step 6 is vulnerable to numerical instability, specifically the sensitivity of the linear system (4) to the condition number of the matrix ${\Phi }_t$. The condition number is as follows:

(18)${\kappa }_t={\displaystyle \frac{{\lambda }_H\left({\Phi }_t\right)}{{\lambda }_L\left({\Phi }_t\right)}},$

where ${\lambda }_H\left({\Phi }_t\right)$ and ${\lambda }_L\left({\Phi }_t\right)$ denote, respectively, the highest and lowest absolute values of the eigenvalues of ${\Phi }_t$. A large condition number indicates that the matrix is close to being singular with $\det {\Phi }_t\approx 0$ (Cheney and Kincaid 2012).

Figure 1 shows three example time series ${\nu }_t$ computed using the above procedure, with the elements of vector ${\sigma }_{\pi t}$ computed using (14). The three time series use log returns from $N=28$ stocks selected from components of the Russell 3000 (RS28), the STOXX Europe 600 (EU28), and the Dow Jones (DJ28) indices. The selection of stocks is discussed in Section 3.1. A historical window of $M=2500$ days was used. The time series illustrates the following issues:

• Ia.. The presence of “spikes”, indicating periods when the matrix ${\Phi }_t$ undergoes rapid changes in the condition number over time.

• Ib.. The investigation of the time series may begin during a period of the rapidly changing condition number (as in the RS28 time series).

• Ic.. During “normal” periods (when the daily change in the condition number ${\Phi }_t$ is random and “small”), the value of ${\nu }_t$ has random behavior.

• Id.. There may be long-term trends in ${\nu }_t$ (seen in all three, but most notably DJ28).

• Ie.. There may be discontinuous changes in the “baseline” value of ${\nu }_t$ associated with rapid changes in the condition number (as in the DJ28 time series).

As illustrated in Figure 2, these spikes always occurred over a multi-day period—the prominent spike in the DJ28 time series occurred from 25 September to 1 October 2018. Thus, regularization techniques that attempt to set ${\nu }_t={\nu }_{t-1}$—when the condition number of ${\Phi }_t$ is flagged as changing significantly—will generally not work.

2.4. Regularization of the Matrix ${\Phi }_t$

Developing a reliable estimate for ${\nu }_t$ requires (RE.i) a method for identifying when the condition number of the matrix ${\Phi }_t$ changes significantly and (RE.ii) implementing a procedure to reduce the change appropriately at those times. The estimation procedure must preserve behaviors Ic and Id while appropriately addressing Ia, Ib, and Ie. Identifying when the condition number ${\kappa }_t$ of ${\Phi }_t$ changes significantly requires a characterization of “what value of ${\kappa }_t-{\kappa }_{t-1}$ is too large?”.

There are a variety of smoothing methods possible for addressing point RE.ii that indirectly control the rate of change of the condition number. One approach is to solve (4) as a least-squares minimization problem with penalty (regularization) terms, i.e.,

(19)${\overline{x}}_t=\underset{z}{arg\; min}\left(\left|\right|{\Phi }_{t}z-{\mu }_t{\left|\right|}_2^2+\left|\right|{\Gamma }_{t}z{\left|\right|}_2^2\right),$

for some choice of the matrix ${\Gamma }_t$, which defines the impact of the regularization on the solution.² An example of such regularization is the minimization, as follows:

(20)${\overline{x}}_t=\underset{z_1,\dots ,z_N}{arg\; min}\left[\left|\right|{\Phi }_{t}z-{\mu }_t{\left|\right|}_2^2+{\gamma }_{1t}{(z_1+{\mu }_{\pi (t-1)})}^2+{\gamma }_{2t}\left(\sum _{k=2}^N{z}_k^2\right)\right],$

where ${\mu }_{\pi (t-1)}$ is the solution to the SRR computed for the previous day and the coefficients ${\gamma }_{1t},{\gamma }_{2t}$ are used to weigh the regularization components. The first penalty term controls the rate of change of the SRR; the second penalty term penalizes solutions with large variance (7). The penalty terms must be time-dependent to avoid regularization when ${\kappa }_t$ does not change significantly. We investigated approach (20). Modeling ${\gamma }_{1t}$ and ${\gamma }_{2t}$ led us to a four-parameter model whose parameters required adjustments with the choice of the asset universe.

We propose a method with several advantages. This method addresses both points RE.i and RE.ii simultaneously and provides insight into the factors controlling ${\kappa }_t$. It requires only a single regularization parameter and is computationally efficient. We consider the singular value decomposition of the matrix ${\Phi }_t$, as follows:

(21)${\Phi }_t=U_t{D}_t{V}_t^{\mathrm{T}}$

where $U_t$ and $V_t$ are orthogonal matrices in $R^{N\times N}$, and $D_t=\mathrm{diag}(d_1^{\left(t\right)},\dots ,d_n^{\left(t\right)})$ consists of the singular values ordered in such a way that $d_1^{\left(t\right)}>d_2^{\left(t\right)}...>d_n^{\left(t\right)}$. The solution to (4) is equivalent to solving the following:

(22a)$\begin{array}{cc}\hfill y& =U_t^{\mathrm{T}}\mu ,\hfill \end{array}$

(22b)$\begin{array}{cc}\hfill z& =D_t^{-1}y,\hfill \end{array}$

(22c)$\begin{array}{cc}\hfill x_t& =V_{t}z.\hfill \end{array}$

Since the condition number of an orthogonal matrix is unity, the ill-conditioning in the solution $x_t$ arises from (22b), specifically, as we show next, due to $z_n=y_n/d_n^{\left(t\right)}$.

Figure 3 plots the time series of singular values $d_j^{\left(t\right)}$, $j=1,\dots ,28$, corresponding to the calculation of ${\nu }_t$ values for EU28 shown in Figure 1. Note the large separations between the group of singular values $d_j^{\left(t\right)}$, $j=2,\dots ,27$, and the largest $d_1^{\left(t\right)}$ and smallest $d_{28}^{\left(t\right)}$ singular values. (This is a consistent finding across all cases studied in Section 3.) The separation of $d_1^{\left(t\right)}$ from the remaining singular values reflects the size difference between the first column of values $11_N$ of the matrix ${\Phi }_t$ and the remaining columns of the submatrix $-{\Sigma }_t$ of ${\Phi }_t$. Figure 3 also plots the time series of (non-zero) singular values $d_j^{(t,\Sigma )}$, $j=1,\dots ,27$ of the submatrix ${\Sigma }_t$. The submatrix ${\Sigma }_t$ is responsible for the intermediate group of singular values $d_j^{\left(t\right)}$, $j=2,\dots ,27$, seen in ${\Phi }_t$.

Variation in the smallest singular value $d_{28}^{\left(t\right)}$ determines the condition number of ${\Phi }_t$ and the spiking behavior shown in ${\nu }_t$. It results when the first column of ${\Phi }_t$, the vector $11_N$, becomes close to the space spanned by the columns of ${\Sigma }_t$. We developed an algorithm to control the time rate of change of the smallest singular value $d_n^{\left(t\right)}$ of ${\Phi }_t$. As the singular values are always positive, we regularize the smallest singular value as follows:

(23)$\begin{array}{cc}\hfill {\overline{d}}_n^{(t,ϵ)}& =\left\{\begin{array}{cc}min\left(d_n^{\left(t\right)},(1+ϵ){\overline{d}}_n^{(t-1,ϵ)}\right)\hfill & \mathrm{if}{d}_n^{\left(t\right)}\ge {\overline{d}}_n^{(t-1,ϵ)},\hfill \\ max\left(d_n^{\left(t\right)},(1-ϵ){\overline{d}}_n^{(t-1,ϵ)}\right)\hfill & \mathrm{if}{d}_n^{\left(t\right)}<{\overline{d}}_n^{(t-1,ϵ)},\hfill \end{array}\right.\hfill \\ \hfill {\overline{d}}_n^{(0,ϵ)}& =d_n^{\left(0\right)}.\hfill \end{array}$

Figure 4 illustrates the sensitivity of the regularization of $d_n^{\left(t\right)}$ to the value $ϵ$; too small a value of $ϵ$ leads to too much regularization, too large a value does not adequately smooth the spikes in $d_n^{\left(t\right)}$.

Let ${\overline{D}}_t^{\left(ϵ\right)}$ denote the matrix D in (21), with the singular value $d_n^{\left(t\right)}$ replaced by the regularized value ${\overline{d}}_n^{(t;ϵ)}$. Let ${\overline{\kappa }}_t^{\left(ϵ\right)}$ denote the condition number of the regularized matrix ${\overline{\Phi }}_t^{\left(ϵ\right)}\equiv U_t{\overline{D}}_t^{\left(ϵ\right)}{V}_t^{\mathrm{T}}$. Figure 4 compares ${\kappa }_t$ against ${\overline{\kappa }}_t^{\left(ϵ\right)}$ for $ϵ\in \{0.001,0.005,0.009\}$. Of course, the behavior of ${\overline{\kappa }}_t^{\left(ϵ\right)}$ relative to ${\kappa }_t$ must reflect the behavior of ${\overline{d}}_n^{(t,ϵ)}$ relative to $d_n^{\left(t\right)}$.

Let ${\overline{\nu }}_t^{\left(ϵ\right)}$ denote the solution obtained from (22a)–(22c) using ${\overline{D}}_t^{\left(ϵ\right)}$ in (22b). Figure 5 compares the results of ${\nu }_t$ against ${\overline{\nu }}_t^{\left(ϵ\right)}$ for $ϵ\in \{0.001,0.005,0.009\}$. The results for the selected EU28 stocks indicate that, of the three, the choice $ϵ=0.009$ provides the best regularization. The plot also indicates the ease with which a choice for $ϵ$ can be made for a particular selection of stocks. The figure also shows details of the behavior of the regularized values ${\overline{\nu }}_t^{\left(ϵ\right)}$ from 1 February 2020 to 1 July 2020. The regularized values “respond” to the extreme volatility expressed by ${\nu }_t$, but in a highly muted manner. There is no visible difference between the $ϵ=0.005$ and $ϵ=0.009$ regularizations over this time period.

The EU28 dataset illustrates two features reflecting the fact that the condition number of ${\Phi }_t$ can change its “baseline” value over time:

1. As the EU28 time series is first sampled at the start of 2011 during a period over which ${\kappa }_t$ is rapidly changing, and the baseline of ${\nu }_t$ appears to be changing, the strength of the regularization determines how soon the regularized values can “settle” to the appropriate condition number and new baseline. In contrast, if the time series is first sampled during a period in which ${\kappa }_t$ is not changing significantly, this initial “transient” behavior seen in the EU28 dataset would not appear.

2. Beginning in 2018, ${\kappa }_t$ begins to change, with $d_{28}^{\left(t\right)}$ decreasing from a baseline value of roughly 0.0025 to a new baseline of roughly 0.001 at the start of 2021. The regularization must accommodate such a baseline change without “looking into the future”, but strongly correct for extreme spiking behavior, as seen in 2020.

Although the value $ϵ=0.009$ produces the best results in Figure 4, it produces only a slight improvement over the choice $ϵ=0.005$. For the other asset groups discussed in Section 3, $ϵ=0.005$ produces the overall best results. For consistency, we proceed with the value $ϵ=0.005$.

2.5. Secondary Regularization of the SRR

The regularization in (23) can be directly applied to the time series ${\overline{\nu }}_t^{\left(ϵ\right)}$, resulting in a further smoothed SRR value, as follows:

(24)$\begin{array}{cc}\hfill {\widehat{\nu }}_t^{(ϵ,{\delta }_{\nu })}& =\left\{\begin{array}{cc}min\left({\overline{\nu }}_t^{\left(ϵ\right)},(1+{\delta }_{\nu }){\widehat{\nu }}_{t-1}^{(ϵ,{\delta }_{\nu })}\right)\hfill & \mathrm{if}{\overline{\nu }}_t^{\left(ϵ\right)}\ge {\widehat{\nu }}_{t-1}^{(ϵ,{\delta }_{\nu })},\hfill \\ max\left({\overline{\nu }}_t^{\left(ϵ\right)},(1-{\delta }_{\nu }){\widehat{\nu }}_{t-1}^{(ϵ,{\delta }_{\nu })}\right)\hfill & \mathrm{if}{\overline{\nu }}_t^{\left(ϵ\right)}<{\widehat{\nu }}_{t-1}^{(ϵ,{\delta }_{\nu })},\hfill \end{array}\right.\hfill \\ \hfill {\widehat{\nu }}_0^{(ϵ,{\delta }_{\nu })}& ={\overline{\nu }}_0^{\left(ϵ\right)}.\hfill \end{array}$

Figure 6 shows the effect of this secondary smoothing using the values ${\delta }_{\nu }\in \{{10}^{-4},{10}^{-5}\}$. As ${\delta }_{\nu }$ decreases further, the time series ${\widehat{\nu }}_t^{(ϵ,{\delta }_{\nu })}$ becomes increasingly linear, losing time-dependent details.

We proceed using the secondary smoothing parameter value ${\delta }_{\nu }={10}^{-5}$.

2.6. Comparison of Solution Methods (14) and (15)

Equations (14) and (15) present alternate methods for approximating the entries ${\sigma }_{jk}$ of the matrix $\Sigma$. The results discussed above have used (14). Here, we briefly illustrate the differences observed between the two methods using the EU28 dataset.

Figure 7 compares the EU28 time series ${\nu }_t$ and ${\widehat{\nu }}_t^{(ϵ,{\delta }_{\nu })}$ for ${\sigma }_{jk}$ obtained via (14) and (15). Using (15) produces an unregularized matrix ${\Phi }_t$ that generally has significant changes in ${\kappa }_t$ more frequently than that obtained using (14). As a result, we apply the regularization (23) to each singular value $d_k^{\left(t\right)}$, $k=1,\dots ,N$, producing the regularized matrix ${\overline{D}}_t^{\left(ϵ\right)}$ and corresponding regularized solution ${\overline{\nu }}_t^{\left(ϵ\right)}$. Secondary regularization via (24) produces the solution ${\widehat{\nu }}_t^{(ϵ,{\delta }_{\nu })}$. The secondary regularized solutions are comparable between the two methods. We have found that, depending on the dataset, either of the two methods may produce secondarily regularized solutions with smaller variations over time.

2.7. Estimation and Regularization of the Total Volatility ${\sigma }_{\pi ,t}$

As noted in Section 2.1, the solution to (4) allows for the determination of the total volatility ${\sigma }_{\pi ,t}$ of the state-price deflator via (7). Let ${\overline{\sigma }}_{\pi k,t}^{\left(ϵ\right)}$, $k=1,\dots ,N-1$, denote the solutions obtained from (22a)–(22c) using ${\overline{D}}_t^{\left(ϵ\right)}$ in (22b). This produces the regularized total volatility ${\overline{\sigma }}_{\pi }^{\left(ϵ\right)}=\sqrt{{\sum }_{k=1}^{N-1}{\left({\overline{\sigma }}_{\pi k,t}^{\left(ϵ\right)}\right)}^2}$. A secondary regularization can also be applied directly to the time series ${\overline{\sigma }}_{\pi k,t}^{\left(ϵ\right)}$, $k=1,\dots ,N-1$, producing the following:

(25)$\begin{array}{cc}\hfill {\widehat{\sigma }}_{\pi k,t}^{(ϵ,{\delta }_{\sigma })}& =\left\{\begin{array}{cc}min\left({\overline{\sigma }}_{\pi k,t}^{\left(ϵ\right)},(1+{\delta }_{\sigma }){\widehat{\sigma }}_{\pi k,t-1}^{(ϵ,{\delta }_{\sigma })}\right)\hfill & \mathrm{if}{\overline{\sigma }}_{\pi k,t}^{\left(ϵ\right)}\ge {\widehat{\sigma }}_{\pi k,t-1}^{(ϵ,{\delta }_{\sigma })},\hfill \\ max\left({\overline{\sigma }}_{\pi k,t}^{\left(ϵ\right)},(1-{\delta }_{\sigma }){\widehat{\sigma }}_{\pi k,t-1}^{(ϵ,{\delta }_{\sigma })}\right)\hfill & \mathrm{if}{\overline{\sigma }}_{\pi k,t}^{\left(ϵ\right)}<{\widehat{\sigma }}_{\pi k,t-1}^{(ϵ,{\delta }_{\sigma })},\hfill \end{array}\right.\hfill \\ \hfill {\widehat{\sigma }}_{\pi k,0}^{(ϵ,{\delta }_{\sigma })}& ={\overline{\sigma }}_{\pi k,0}^{\left(ϵ\right)}.\hfill \end{array}$

The difference in the scale of the values ${\nu }_t$ and ${\sigma }_{\pi k,t}$ necessitates that ${\delta }_{\sigma }>{\delta }_{\nu }$. We found the value ${\delta }_{\sigma }={10}^{-3}$ provided adequate smoothing.

Figure 8 compares the unregularized time series ${\sigma }_{\pi ,t}$ from (7) with the regularized time series ${\widehat{\sigma }}_{\pi ,t}^{(ϵ,{\delta }_{\sigma })}=\sqrt{{\sum }_{k=1}^{N-1}{\left({\widehat{\sigma }}_{\pi k,t}^{(ϵ,{\delta }_{\sigma })}\right)}^2}$. Shown in the figure are comparisons using the solutions (14) and (15) for the matrix $\Sigma$.

3. Results: Empirical Application

A projected use of the SRR is as a discriminator between asset classes, with the implication that an asset class with higher values of ${\nu }_t$ is preferred. We tested the estimated SRR on four different empirical datasets of stocks and included an additional dataset of exchange-traded funds.

3.1. Variation with Asset Type

Each dataset consisted of $N=28$ assets. The first asset group (denoted DJ28) is composed of the stocks in the Dow Jones Industrial Average (excluding the Dow Chemical Company and Visa Inc., for which the historical data (from 20 March 2019 and 18 March 2008, respectively) were too limited).³ The second (denoted SP28), third (denoted RS28), and fourth (denoted EU28) asset groups are composed of subsets of the Standard and Poor’s 500 Index, the Russell 3000 Index, and the Stoxx Europe 600 Index, respectively. The last group (denoted ETF28) consists of a selection of ETFs on US stocks. To avoid selection bias, we utilized a consistent method to select the subset of securities for each of SP28, RS28, EU28, and ETF28. Specifically, we performed the following steps:

• Ordered the assets in the universe in question by their market capitalization;

• Parsimoniously and symmetrically removed the assets with the lowest and highest capitalization until the remaining number was divisible by 28; and

• Picked the assets corresponding to the $(k-0.5)/28$, $k=1,\dots ,28$, percentile assets.

Each dataset consisted of daily asset prices from 3 January 2001 to 29 June 2022. We used an estimation window of ten years (specifically 2500 trading days). Daily SRRs were computed from 3 January 2001 to 29 June 2022. As noted in Section 2.1, we utilized the regularization parameter values $ϵ=0.005$ and ${\delta }_{\nu }={10}^{-5}$.

Figure 9 shows the regularized SRR time series ${\widehat{\nu }}_t^{(ϵ,{\delta }_{\nu })}$ obtained for these five 28-asset groups. Note that changes in the condition numbers of their respective matrices ${\Phi }_t$ over time led to a drift in the “baseline” behavior of ${\nu }_t$ (and, hence, for ${\widehat{\nu }}_t^{(ϵ,{\delta }_{\nu })}$). Figure 10 presents box–whisker summaries of the distribution of values for ${\kappa }_t$, ${\overline{\kappa }}_t^{\left(ϵ\right)}$, ${\nu }_t$, and ${\widehat{\nu }}_t^{(ϵ,{\delta }_{\nu })}$ over this 11.5-year period. The quantile values obtained for ${\widehat{\nu }}_t^{(ϵ,{\delta }_{\nu })}$ are presented in Table 1. The results strongly support the DJ28 group as the preferred investment group, showing positive SRR values over 75% of the 11.5-year period. The EU28 group also demonstrates positive SRR values for over 75% of the 11.5-year period, but its daily SRR return rates exhibit much narrower variation. For the remaining three groups, the regularized SRR values are negative for more than 50% of the 11.5-year time period.

3.2. Dependence on Group Size

We computed the SRR for a group of 1252 assets chosen from the Russell 3000 index. The choice of assets was determined by the requirement that asset prices cover the period from 10 January 2000 to 29 June 2022. Again, daily SRR values were computed from 3 January 2001 to 29 June 2022 using a moving window of 2500 days. Regularization was conducted using the same parameter values as in Section 3.1 with (14) used for the computation of the elements of $\Sigma$.

Figure 11 shows the time series of singular values $d_j^{\left(t\right)}$, $j=1,\dots ,1252$, of the matrix ${\Phi }_t$ for this group. As for the 28-asset groups discussed in Section 3.1, the behavior of the smallest singular value is responsible for the behavior of the condition number of the 2891 × 1252 matrix ${\Phi }_t$. From the behavior of $d_{1252}^{\left(t\right)}$, it is fairly obvious that the matrix ${\Phi }_t$ is both poorly conditioned and has frequent, significant changes in the condition number.

The respective plots in Figure 10 include the distribution of unregularized ${\kappa }_t$ and regularized ${\overline{\kappa }}_t^{\left(ϵ\right)}$ condition numbers seen over this time period for the 1252-asset group. While the regularized condition number is ${\overline{\kappa }}_t^{ϵ}\sim O\left({10}^3\right)-O\left({10}^5\right)$ for the 28-asset groups, the regularized condition number is ${\overline{\kappa }}_t^{ϵ}\sim O\left({10}^7\right)$ for the 1252-asset group. The respective plots in Figure 10 also include the distributions of ${\nu }_t$ and ${\widehat{\nu }}_t^{(ϵ,{\delta }_{\nu })}$ for the RS1252 group. Despite its poorly behaved condition number, the results for RS1252 show improved regularized SRR values compared to the RS28 group. The values of the quartiles of this distribution are also presented in Table 1.

3.3. Behavior of the State-Price Deflator

As noted in Section 2.1, computations of the drift ${\mu }_{\pi ,t}$ and total volatility ${\sigma }_{\pi ,t}$ time series offer critical insight into the trajectory of the state-price deflator ${\pi }_t$. As the SRR ${\nu }_t=-{\mu }_{\pi ,t}$, we consider the behavior of the pair ${\nu }_t$ and ${\sigma }_{\pi ,t}$, using the regularized values ${\widehat{\nu }}_t^{(ϵ,{\delta }_{\nu })}$ and ${\widehat{\sigma }}_{\pi ,t}^{(ϵ,{\delta }_{\sigma })}$.

Comparisons of the distributions of the SRR values ${\widehat{\nu }}_t^{(ϵ,{\delta }_{\nu })}$ for the six asset groups analyzed in Section 3.1 and Section 3.2 are presented in Figure 10. The distributions of the values ${\widehat{\sigma }}_{\pi ,t}^{(ϵ,{\delta }_{\sigma })}$ for the six asset groups are presented in Figure 12. The volatility of the 1252-asset group is significantly larger than that of the 28-asset groups, undoubtedly reflecting the fact that the deflator must now ensure $\mathbb{P}-$martingale behavior for a larger set of correlated assets.

Figure 12 displays the trajectories of the point $\left({\widehat{\sigma }}_{\pi ,t}^{(ϵ,\delta )},{\widehat{\nu }}_t^{(ϵ,{\delta }_{\nu })}\right)$ over this time period for all six asset groups. While there are time periods of linear correlation between ${\widehat{\nu }}_t^{(ϵ,{\delta }_{\nu })}$ and ${\widehat{\sigma }}_{\pi ,t}^{(ϵ,\delta )}$ (with either a positive or negative slope), there are “breakpoints” when the behavior abruptly changes.

As the DJ28 trajectory is the simplest, we explore it in more detail. Figure 13 plots unregularized and regularized time series of ${\nu }_t$ and ${\sigma }_{\pi ,t}$ for this 11.5-year time period. Except for the prominent spike from 24 September to 1 October 2018 (discussed in Section 2.3), the regularized solution smoothed the overall downward trend in ${\nu }_t$, while capturing the prominent “bump” that occurred from 30 October 2015 to 21 September 2018. While the spike resulted in a lowering of the base value of ${\nu }_t$, it correlated with a change from overall downward to flat-to-increasing behavior in ${\sigma }_{\pi ,t}$. There was also a prominent upward “bump” in ${\sigma }_{\pi ,t}$ values from 30 October 2015 to 21 September 2018.

Figure 14 summarizes the trajectories of the unregularized time values $\left({\sigma }_{\pi ,t},{\nu }_t\right)$ and the regularized time values $\left({\widehat{\sigma }}_{\pi ,t}^{(ϵ,{\delta }_{\sigma })},{\widehat{\nu }}_t^{(ϵ,{\delta }_{\nu })}\right)$ for the DJ28 group. The range of axis values is determined by the regularized parameter plot on the right. As a result, seven outlier values from the unregularized parameter plot on the left are not visible. These outlier values correspond to spike events occurring on trading days from 24 September to 2 October 2018. The start and end dates of the spike discontinuity are displayed in the unregularized parameter plot. The regularized parameter plot eliminates the discontinuity and introduces a relaxation period (3 October 2018 to 21 August 2020) to smooth the discontinuity. Outside of this period, the regularized parameter trajectory provides a reasonably accurate smoothing of the unregularized trajectory.

4. Discussion

Our principal concern in this work involved the development of a numerical model to empirically compute the shadow riskless rate described in Section 2.1 for any selected universe of assets. There are two significant challenges to this. The first challenge involves the approximation of the volatilities ${\sigma }_{jk}$ in (1) describing the stochastic effects of the $N-1$ Brownian motions on each risky asset. The method used for estimating these parameters is detailed in Section 2.2. The second challenge involves controlling the ill-conditioned nature of the matrix $\Phi$ of the linear system (4) that uniquely determines the shadow riskless rate. As a rough rule-of-thumb, a condition number of $\kappa (\Phi )={10}^k$ implies a potential loss of k digits of numerical accuracy, in addition to the loss of precision inherent in the LU decomposition (with pivoting) used to solve (4). There is, of course, no way to replace the linear system (4) with a well-conditioned linear system that has the same solution set $\{x^{\left(i\right)},i=1,...\}$ for any given input set $\{{\mu }^{\left(i\right)},i=1,...\}$. Rather, one has to develop an approximate (regularized) solution that controls changes in the condition number of the system. Our method and its advantages are presented in Section 2.4 and Section 2.5.

A side benefit of the empirical work is that it provides insight into the time development of the drift and total volatility of the state-price deflator for the asset group considered. We investigated this in Section 3.3.

Given the relatively brief exploration of datasets in Section 3.1 and Section 3.2, two important avenues require deeper investigation. The first involves investigating the SRR for other asset classes (fixed income, commodities, crypto-assets, foreign currency, etc.) as well as a variety of geopolitical financial markets (China, United Kingdom, Japan, India, Hong Kong, etc.). The second requires a more comprehensive investigation of the scaling of the SRR with the size of the asset group. Of further-reaching importance to the asset-class discriminatory use of the SRR is the question of how well future values (and how far into the future such values) can be predicted. This approach might be accomplished by fitting a historical SRR time series, appropriately modified to be stationary, to an ARMA-GARCH model and developing multi-day, out-of-sample forecasts. Value-at-risk backtesting can then be applied to the out-of-sample projections to determine the loss of predictive accuracy with increasing forecast time. Since the assumption of asset prices governed by geometric Brownian motion fundamentally conflicts with observed price return distributions, as noted in the theoretical work by Rachev et al. (2017), it will ultimately be necessary to develop a procedure to compute SRR values for risky assets following more realistic price processes, such as those described by (i) jump diffusion, (ii) diffusion with local volatility, and (iii) geometric fractional Brownian or Rosenblatt motion.

Prior to pursuing these investigations, our work suggests a second avenue for defining a shadow riskless rate, whose numerical solution would likely be much better conditioned. We start by noting that the column eigenvectors $w_{\circ j}$, $j=1,\dots ,N$, obtained from the PCA analysis of the historical return data for the N underlying risky assets, are normalized. Column vector $w_{\circ j}$ can be interpreted as a unique (and orthogonal) set of N positions. Applying these positions to the N underlying risky assets produces a composite asset ${\mathcal{R}}^{\left(j\right)}$ with a return $r^{\left(j\right)}={\sum }_{i=1}^N{w}_{ij}{r}_i$, $j=1,\dots ,N$, where $r_i$ represents a historical return for underlying asset i. Assume N is large and let $1<K\ll N$.

Consider the following algorithm used to estimate a low-variance rate (Algorithm 1).

Algorithm 1: SRR based upon PCA and mean-variance portfolio optimization

Find the minimum mean-variance portfolio ${\sigma }_p^{\left(K\right)},r_p^{\left(K\right)}$ for the composite assets           ${\mathcal{R}}^{\left(1\right)},\dots ,{\mathcal{R}}^{\left(K\right)}$ $j=K$ continue = true while (continue)       $j=j+1$       Find the minimum mean-variance portfolio ${\sigma }_p^{\left(j\right)},r_p^{\left(j\right)}$ for the composite assets           ${\mathcal{R}}^{\left(1\right)},\dots ,{\mathcal{R}}^{\left(j\right)}$       if ((${\sigma }_p^{\left(j\right)}-{\sigma }_p^{(j-1)}>{\mathrm{tol}}_{\sigma }$) or ($r_p^{\left(j\right)}-r_p^{(j-1)}>{\mathrm{tol}}_r$)) continue = false end $r=r_p^{\left(j\right)}$ ${\sigma }_r={\sigma }_p^{\left(j\right)}$

This final minimum mean-variance portfolio provides a set of weights $q_r=\{q_1^{\left(j\right)},\dots ,q_j^{\left(j\right)}\}$ satisfying ${\sum }_{i=1}^j{q}_i^{\left(j\right)}=1$.⁴

As the eigenvectors $w_{\circ j}$ are ordered according to the amount of variance they explain in the data, gradually adding composite assets ${\mathcal{R}}^{\left(j\right)}$, $j=K,K+1,...$ should slowly decrease both ${\sigma }_p^{\left(j\right)}$ and $r_p^{\left(j\right)}$ until a value of j is reached, such that successive additions merely add “noise” that will increase either ${\sigma }_p^{\left(j\right)}$ or $r_p^{\left(j\right)}$. The return value r from Algorithm 1 is the minimum return possible (under the variance risk measure) for this asset class, with a minimum standard deviation ${\sigma }_r$. We propose that r can provide an empirical proxy to a minimum-risk rate for the class of assets under consideration. The weights $q_r$ applied to the composite assets ${\mathcal{R}}^{\left(1\right)},\dots ,{\mathcal{R}}^{\left(j\right)}$ define the minimum-variance asset whose drift component is given by r.

In contrast, let $\left(r^{\left(N\right)},{\sigma }_r^{\left(N\right)}\right)$ denote the minimum mean-variance portfolio achieved by considering the full set N of underlying assets, each having historical return $r_i$. The value of ${\sigma }_r$ should be smaller than that of ${\sigma }_r^{\left(N\right)}$, with the question being, “can it be significantly smaller?”.

Author Contributions

Conceptualization, S.T.R.; methodology, W.B.L., D.L., J.P. and S.T.R.; software, D.L., W.B.L., J.P. and Y.H.; validation, W.B.L., D.L. and S.T.R.; investigation, W.B.L., D.L., J.P. and Y.H.; data curation, D.L. and J.P.; writing—original draft preparation, D.L. and J.P.; writing—review and editing, W.B.L. and F.J.F.; visualization, D.L. and W.B.L.; supervision, S.T.R. and F.J.F. All authors have read and agreed to the published version of the manuscript.

Data Availability Statement

The asset data used in this study were provided under a license from Bloomberg Professional Services. A copy of the SRR computational code is available upon request from the corresponding author.

Conflicts of Interest

Author Jiho Park was employed by the company Citigroup. The remaining authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest. The [company Citigroup in affiliation] had no role in the design of the study; in the collection, analyses, or interpretation of data; in the writing of the manuscript, or in the decision to publish the results.

Footnotes

Notes

1. Attempts to use robust linear regression with M-estimation produced a $\Sigma$ matrix whose condition number and frequency of change of condition number were even worse.

2. The solution to (19) is equivalent to finding the solution $x_t$ that minimizes the $L_2$-norm of the error vector ${ϵ}_t$ of the regression problem ${\mu }_t={\Phi }_t{x}_t+{ϵ}_t$. Approach (19) is known as Tikhonov regularization (Tikhonov and Glasko 1965) when applied to integral equations, and known as ridge regression (Hoerl and Kennard 1970) when applied to finite-dimensional regression problems.

3. (Hence, the rationale for setting $N=28$).

4. We assume long-only mean-variance optimization.

5. The matrix $X^{\mathrm{T}}X$ is positive semidefinite.

Footnotes

Disclaimer/Publisher’s Note: The statements, opinions and data contained in all publications are solely those of the individual author(s) and contributor(s) and not of MDPI and/or the editor(s). MDPI and/or the editor(s) disclaim responsibility for any injury to people or property resulting from any ideas, methods, instructions or products referred to in the content.

Figures and Table

Figure 1. Plots of the SRR [Forumla omitted. See PDF.] computed using steps 1 to 6, with (14) used in step 4. Note that the large positive spike in the plotted time series for DJ28 is truncated to 11% of its true height.

Figure 2. Plots of the SRR [Forumla omitted. See PDF.] computed using (5).

Figure 3. Time series of (top) the singular values [Forumla omitted. See PDF.], [Forumla omitted. See PDF.] of the matrix [Forumla omitted. See PDF.] and (bottom) the singular values [Forumla omitted. See PDF.], [Forumla omitted. See PDF.] of the matrix [Forumla omitted. See PDF.] for the RS28 dataset. Each singular value, ordered from largest to smallest, is represented by a colored line.

Figure 4. The EU28 time series (top) [Forumla omitted. See PDF.] and [Forumla omitted. See PDF.]; (bottom) [Forumla omitted. See PDF.] and [Forumla omitted. See PDF.] for the regularization values [Forumla omitted. See PDF.].

Figure 4. The EU28 time series (top) [Forumla omitted. See PDF.] and [Forumla omitted. See PDF.]; (bottom) [Forumla omitted. See PDF.] and [Forumla omitted. See PDF.] for the regularization values [Forumla omitted. See PDF.].

Figure 5. (top)The EU28 times series [Forumla omitted. See PDF.] and [Forumla omitted. See PDF.]; (bottom) a detailed view, from 1 February 2020 to 1 July 2020, of these two series for the regularization values [Forumla omitted. See PDF.].

Figure 6. The EU28 time series [Forumla omitted. See PDF.], [Forumla omitted. See PDF.] and [Forumla omitted. See PDF.] for the regularization values [Forumla omitted. See PDF.].

Figure 7. The EU28 time series [Forumla omitted. See PDF.] and [Forumla omitted. See PDF.] for [Forumla omitted. See PDF.] obtained via (top) (14) and (bottom) (15).

Figure 8. The EU28 time series [Forumla omitted. See PDF.] and [Forumla omitted. See PDF.] for [Forumla omitted. See PDF.] obtained via (top) (14) and (bottom) (15).

Figure 8. The EU28 time series [Forumla omitted. See PDF.] and [Forumla omitted. See PDF.] for [Forumla omitted. See PDF.] obtained via (top) (14) and (bottom) (15).

Figure 9. The time series [Forumla omitted. See PDF.] for the five 28-asset groups.

Figure 10. Box–whisker summaries of the time-series distributions of condition numbers (top-left) [Forumla omitted. See PDF.], (top-right) [Forumla omitted. See PDF.], and SRR values (bottom-left) [Forumla omitted. See PDF.], (bottom-right) [Forumla omitted. See PDF.], for the five 28-asset groups and the RS1252 group. The summaries for SRR are plotted using the y-axes limits determined by the [Forumla omitted. See PDF.] box–whisker plots; hence, many outliers for the [Forumla omitted. See PDF.] box-whisker plots are not shown.

Figure 11. The time series of the singular values [Forumla omitted. See PDF.], [Forumla omitted. See PDF.] of the matrix [Forumla omitted. See PDF.] for the 1252-asset group.

Figure 12. (left) Box–whisker summaries of the regularized time-series distributions of SS values for the five 28-asset groups and the RS1252 group. (right) The trajectories of the point [Forumla omitted. See PDF.] for the 1252-asset and 28-asset groups.

Figure 13. Time series for the regularized and unregularized values of (top) [Forumla omitted. See PDF.] and (bottom) [Forumla omitted. See PDF.] for the DJ28-asset group.

Figure 14. Trajectory of (left) [Forumla omitted. See PDF.] and (right) [Forumla omitted. See PDF.] for the DJ28 asset group, with significant time points indicated.

Table 1

Quantile values of ${\widehat{\nu }}_t^{(ϵ,{\delta }_{\nu })}$ from Figure 10.

Quantile	DJ28	ETF27	EU28	RS28	SP28	RS1252
$P_{75}$	4.9 $\times {10}^{-3}$	$1.4\times {10}^{-4}$	$6.3\times {10}^{-4}$	$-2.5\times {10}^{-4}$	$9.3\times {10}^{-5}$	$1.1\times {10}^{-3}$
$P_{50}$	$3.7\times {10}^{-3}$	$-3.8\times {10}^{-4}$	$1.9\times {10}^{-4}$	$-9.5\times {10}^{-4}$	$-2.0\times {10}^{-4}$	$4.9\times {10}^{-4}$
$P_{25}$	$-7.1\times {10}^{-5}$	$-9.4\times {10}^{-4}$	$4.4\times {10}^{-5}$	$-2.2\times {10}^{-3}$	$-1.5\times {10}^{-3}$	$6.2\times {10}^{-5}$

Appendix A. Principal Component Analysis

Consider the data matrix $X\in {\mathbb{R}}^{M\times N}$, where each row vector $x_{m\circ }$ of X represents the values obtained from a sample m, where $m=1,\dots ,M$, with each sample consisting of N elements. We assume each column vector $x_{\circ j}$ of X satisfies ${\displaystyle \frac{1}{M}}{\sum }_{m=1}^M{x}_{mj}=0$ (i.e., the sample mean of each column has been shifted to zero). Let column $w_{\circ j}$, $j=1,\dots ,N$, of the matrix $W\in {\mathbb{R}}^{N\times N}$ denote the orthonormal eigenvectors of $X^{\mathrm{T}}X$. Then, $X^{\mathrm{T}}XW=W\Lambda$, where $\Lambda =\mathrm{diag}({\lambda }_1,\dots ,{\lambda }_N)$ contains the eigenvalues. The eigenvalue–eigenvector pairs are ordered such that ${\lambda }_1\ge {\lambda }_2\ge ...\ge {\lambda }_N\ge 0$.5 The matrix $P\in {\mathbb{R}}^{M\times N}$ satisfying $P=XW$ has the property that each element $p_{mj}$ represents the projection of the row vector $x_{m\circ }$ onto the eigenvector $w_{\circ j}$. Thus, the components $p_{mj}$, $m=1,\dots ,M$, of the column vector $p_{\circ j}$ of P, satisfying $p_{\circ j}=X{w}_{\circ j}$, represent the projections of the vectors $x_{m\circ }$ (i.e., the results of experiment m) on the eigenvector $w_{\circ j}$. The vectors $p_{\circ j}=X{w}_{\circ j}$ are known as the principal components of the analysis. The sample covariance between two different principal components is zero due to the orthogonality of the eigenvectors:(A1)$$\mathrm{Cov}(p_{\circ i},p_{\circ j})={\left(X{w}_{\circ i}\right)}^{\mathrm{T}}\left(X{w}_{\circ j}\right)={\lambda }_j{w}_{\circ i}^{\mathrm{T}}{w}_{\circ j}={\lambda }_j{\delta }_{ij}.$$

Thus, the principal component covariance matrix $C_P$ is $P^{\mathrm{T}}P=\Lambda$. The principal component analysis transforms data to coordinates $w_{\circ j},j=1,\dots ,N$, which diagonalize the empirical covariance matrix.

Given that $X^{\mathrm{T}}XW=W\Lambda =W\Lambda W^{\mathrm{T}}W$, and considering that the orthonormal columns $w_{\circ j}$ of W span ${\mathbb{R}}^N$, we can express the covariance matrix $C_X$ of the empirical data as follows:(A2)$$C_X=X^{\mathrm{T}}X=W\Lambda W^{\mathrm{T}}=\left(W\sqrt{\Lambda }\right){\left(W\sqrt{\Lambda }\right)}^{\mathrm{T}}.$$

Thus, for the case where $M>N$, the principal component analysis identifies a reduced space ${\mathbb{R}}^{N\times N}$ spanned by an orthonormal set of vectors $w_{\circ j}$, $j=1,\dots ,N$. The spanning vectors $w_{\circ j}$ are eigenvectors of the data matrix $X^{\mathrm{T}}X$, ordered such that the first principal component $p_{\circ 1}=X{w}_{\circ 1}$ explains the greatest variance of the data; the second principal component $p_{\circ 2}=X{w}_{\circ 2}$ explains the next largest amount of variance, etc.

References

Altinoglu, Levent. A Theory of Safe Asset Creation, Systematic Risk, and Aggregate Demand; Finance and Economics Discussion Series 2023-062 Board of Governors of the Federal Reserve System: Washington, DC, 2023.

Aoki, Kosuke; Nakajima, Tomoyuki; Nikolov, Kalin. Safe asset shortages and asset price bubbles. Journal of Mathematical Economics; 2014; 3, pp. 164-74. [DOI: https://dx.doi.org/10.1016/j.jmateco.2014.05.005]

Black, Fischer. Capital market equilibrium with restricted borrowing. The Journal of Business; 1972; 45, pp. 444-55. [DOI: https://dx.doi.org/10.1086/295472]

Bocola, Luigi; Lorenzoni, Guido. Risk-sharing externalities. Journal of Political Economy; 2023; 131, pp. 595-632. [DOI: https://dx.doi.org/10.1086/722088]

Caballero, Ricardo J.; Simsek, Alp. A risk-centric model of demand recessions and speculation. The Quarterly Journal of Economics; 2020; 135, pp. 1493-566. [DOI: https://dx.doi.org/10.1093/qje/qjaa008]

Caballero, Ricardo J.; Simsek, Alp. A model of endogenous risk intolerance and LSAPs: Asset prices and aggregate demand in a COVID-19 shock. The Review of Financial Studies; 2021; 34, pp. 5522-80. [DOI: https://dx.doi.org/10.1093/rfs/hhab036]

Caballero, Ricardo J.; Farhi, Emmanuel. A model of the safe asset mechanism (SAM): Safety traps and economic policy. Working Paper 18737; National Bureau of Economic Research: Cambridge, MA, 2013.

Caballero, Ricardo J.; Farhi, Emmanuel. The safety trap. The Review of Economic Studies; 2018; 85, pp. 223-74.

Caramp, Nicolas. Sowing the seeds of financial crises: Endogenous asset creation and adverse selection. The Review of Economic Studies; 2024; [DOI: https://dx.doi.org/10.1093/restud/rdae092]

Cheney, E. Ward; Kincaid, David R. Numerical Mathematics and Computing; Cengage: Boston, 2012.

Duffie, Darrell. Dynamic Asset Pricing Theory; 3rd ed. Princeton University Press: Princeton, 2001.

Eggertsson, Gauti B.; Krugman, Paul. Debt, deleveraging, and the liquidity trap: A Fisher–Minsky–Koo approach. The Quarterly Journal of Economics; 2012; 127, pp. 1469-513. [DOI: https://dx.doi.org/10.1093/qje/qjs023]

Gorton, Gary. The history and economics of safe assets. Annual Review of Economics; 2017; 9, pp. 547-86. [DOI: https://dx.doi.org/10.1146/annurev-economics-033017-125810]

Gorton, Gary; Lewellen, Stefan; Metrick, Andrew. The safe-asset share. American Economic Review; 2012; 102, pp. 101-6. [DOI: https://dx.doi.org/10.1257/aer.102.3.101]

Gorton, Gary; Laarits, Toomas; Muir, Tyler. Mobile collateral versus immobile collateral. Journal of Money, Credit and Banking; 2022; 54, pp. 1673-703. [DOI: https://dx.doi.org/10.1111/jmcb.12917]

Gourinchas, Pierre-Olivier; Jeanne, Olivier. Global safe sssets. BIS Working Papers 399; Bank for International Settlements: Basel, 2012.

Hoerl, Arthur E.; Kennard, Robert W. Ridge regression: Biased estimation for nonorthogonal problems. Technometrics; 1970; 12, pp. 55-67. [DOI: https://dx.doi.org/10.1080/00401706.1970.10488634]

Holmstrom, Bengt. Understanding the role of debt in the financial system. BIS Working Papers 479; Bank for International Settlements: Basel, 2015.

Rachev, Svetlozar T.; Stoyanov, Stoyan V.; Fabozzi, Frank J. Financial markets with no riskless (safe) asset. International Journal of Theoretical and Applied Finance; 2017; 20, 1750054. [DOI: https://dx.doi.org/10.1142/S0219024917500546]

Tikhonov, Andrei Nikolaevich; Glasko, Vladlen Borisovich. Use of the regularization method in non-linear problems. USSR Computational Mathematics and Mathematical Physics; 1965; 5, pp. 93-107. [DOI: https://dx.doi.org/10.1016/0041-5553(65)90150-3]