Empirical estimation of realized higher-order (joint) cumulants and their practical issues

As stock returns are not fully explained by the classical asset pricing model, various models have been developed to explain anomalies. One approach extends the mean-variance framework by incorporating higher-order moments of the return distribution, such as skewness and kurtosis. Kraus and Litzenberger (1976), Kalev et al. (2019), and Harvey and Siddique (2000) develop an asset pricing model incorporating coskewness between returns and an index. Brunnermeier et al. (2007), Mitton and Vorkink (2007), and Barberis and Huang (2008) suggest that skewness or idiosyncratic skewness of returns affects the expected returns of assets. Dittmar (2002), Cvitanić et al. (2008), Chabi-Yo (2012), and Duan and Zhang (2014) suggest that kurtosis also affects the expected returns. These theoretical suggestions are supported by extensive empirical studies that provide evidence of the relationship between the non-normality of the return distribution and the excess returns of financial assets (Friesen et al., 2012; Amaya et al., 2015; Stilger et al., 2017; Langlois, 2020; Fan et al., 2022).

In empirical studies, higher order moments are estimated through various methods. In particular, Amaya et al. (2015), Ahadzie and Jeyasreedharan (2020), and Dai et al. (2023) define realized k-th moments as the sum of k-th powers of sub-period returns, $∑i=1N(Sti−Sti−1)k,0=t0<t1<⋯<tN=T$ with $k≥2$⁠, where $Sti$ is the price of a security at time $ti$⁠. However, the above estimators (hereafter, conventional realized moments) are biased estimators of the k-th true moments $E0[(ST−S0)k]$ for $k≥3$⁠. To overcome this weakness, Neuberger (2012), Bae and Lee (2021), Fukasawa and Matsushita (2021), and Bae (2022) propose realized higher-order (joint) cumulants and show that these provide unbiased estimators of cumulants for a specific period ^[1]. The estimators are developed based on an extended information set containing variance and other (co)moments of returns up to time T. As an example, the third moment of a period consists not only of the third power of sub-period returns but also of the relation between variance and sub-period returns. Similarly, the fourth cumulant of a period is affected by the volatility structure and the relation between skewness and returns, in addition to the fourth power of the sub-period returns.

In this study, we discuss practical issues in estimating realized cumulants and their applications. First, we provide examples to obtain (co)moments using option prices. We can use them as implied (co)moments themselves or as extended information to obtain practical realized (joint) cumulants. Second, through simulation analysis, we confirm that implied (joint) cumulants are highly correlated with the true ex ante (joint) cumulants, although the average values of implied cumulants differ from those of the true cumulants. We also find that practically obtained realized cumulants are highly correlated with the true realized (joint) cumulants. Lastly, we estimate realized (co)skewness and kurtosis using the S&P 500 index, individual stock prices, and their options, and examine the predictability. As a result, we find that the realized (co)skewness and kurtosis are predicted by implied and lagged realized ones. According to Jiang and Tian (2005), implied variance has information content about realized variance ^[2]. Similarly, we can expect that implied skewness and kurtosis also have predictability for realized skewness and kurtosis, and the above empirical results confirm this ^[3]. Considering the efficiency of the options market, these results suggest that realized moments appropriately represent the characteristics of realized return over a period.

Prior studies, including Kraus and Litzenberger (1976), Kalev et al. (2019), and Harvey and Siddique (2000), argue that investor utility depends not only on the mean and variance of returns but also on their skewness and kurtosis. According to these studies, optimal investment decisions require information on the future skewness and kurtosis of returns. Bakshi et al. (2003) propose methods to extract implied skewness and kurtosis from option prices; however, these quantities are obtained under the risk-neutral (Q) measure, and thus may differ from the moments under the physical (P) measure, which are relevant to investors' utility. In this study, we show that Q-measure skewness and kurtosis are positively correlated with their P-measure counterparts, which implies that skewness and kurtosis extracted from options can still be used meaningfully in investment decision-making.

The rest of the paper is organized as follows. Section 2 outlines the estimation of realized (joint) cumulants based on Bae and Lee (2021). Section 3 examines the validity of the practically obtained realized (joint) cumulants through simulation analysis, while Section 4 investigates their predictors through empirical analysis using the S&P 500 index, individual stock prices, and their option prices.

Neuberger (2012) and Bae and Lee (2021) propose Eqs. (6)–(10) for realized third- and fourth-order (joint) cumulants, which satisfy Eqs. (1)–(5) for the martingale processes $S1,t$ and $S2,t$⁠.

with

where $ΔXtj=Xtj−Xtj−1$ and $ΔXi,tj=Xi,tj−Xi,tj−1$ for any process $X$ and asset $i$⁠, given a partition ${t0,t1,…,tN}$ on $[0,T]$ with $0=t0≤t1≤t2≤…≤tN=T$⁠.

Additionally, the authors show that $∑j=1N(ΔS1,tj)k1(ΔS2,tj)k2$ for $k1+k2≥3$ are biased estimators of the true comoments, $E0[(S1,T−S1,0)k1(S2,T−S2,0)k2]$⁠. Specifically, according to Eq. (8), conventional third comoment $∑j=1N(ΔS1,tj)2(ΔS2,tj)$ is a biased estimator of the true third comoment. To obtain unbiased estimates, additional terms $(ΔS1,tj)(ΔCtj)$ and $(ΔS2,tj)(ΔV1,tj)$ are required, which capture the contagion effect and the interdependence between different assets. Similarly, the fourth cumulant requires the terms $(ΔS1,j)2(ΔV1,tj),(ΔV1,tj)2$⁠, and $ΔS1,j(ΔW1,tj)$⁠, which capture the volatility clustering effect and the correlation between skewness and returns.

However, there are several practical issues in estimating the realized (joint) cumulants. First, since the equations above are derived under the assumption that the price of each security is a martingale, we can use forward prices for each $Si,j$⁠. Without loss of generality, we assume $Si,0$ equals one. In other words, hereafter, $Si,j$ denotes $Si,j/Si,0$⁠, so that $Si,T−Si,0$ is the arithmetic return between time 0 and $T$⁠.

Second, to estimate the realized fourth cumulant, we need the second and third moments. Although we cannot directly observe these moments, we can obtain implied moments of $Si,T−Si,j$ by adopting the methods of Bakshi and Madan (2000) and Bakshi et al. (2003) as follows:

where $Pj(x)$ is the forward price of a European put option at time $j$ with exercise price $x$ and maturity $T$⁠. Similarly, $Cj(x)$ is the forward price of a European call option.

Third, to estimate the realized joint cumulant, we need the comoments between two assets. The comoments are neither directly observable nor obtainable from individual European options alone. Instead, we can obtain the implied comoments using certain exotic options. For example, if there exists one basket option on $a1S1,T+a2S2,T$ and two individual stock options on $S1,T$ and $S2,T$⁠, we can obtain information about $varj(a1S1,T+a2S2,T)$⁠, $varj(S1,T)$⁠, and $varj(S2,T)$ using Eq. (16), and then the covariance of the two assets can be derived as follows:

If we have additional basket options on $a3S1,T+a4S2,T$⁠, where $a3$ and $a4$ are non-zero constants with $a1a4≠a2a3$⁠, we can obtain implied third moments for $S1,T$⁠, $S2,T$⁠, $a1S1,T+a2S2,T$⁠, and $a3S1,T+a4S2,T$⁠. Using this information, we can obtain the implied third cumulant between the two assets as follows:

where $TMjimp$ represents the implied third moment at time $j$⁠.

Even without the exotic options described above, we can obtain the covariance between an index and its component by adopting the method of Kempf et al. (2015) ^[4], provided that index options are available. Since an index is a portfolio of individual stocks, it is possible to extract joint cumulants between the index and an individual asset under certain assumptions about the relationships between the asset and the index. One of the assumptions is that asset returns follow an index model with time-varying $α$ and $β$⁠. Under this assumption, at each time $tj$⁠, the relation between the stock index (⁠$SM,T$⁠) and the individual stock price (⁠$Si,T$⁠) can be represented as follows:

The second assumption is that the ratio of systematic risk to total risk is the same for all securities at each time $tj$ denoted by $ρj$⁠. With these two assumptions, we obtain $βi,j2VM,j=ρjVi,j$ or

where $Vi,j=varj(Si,T)$ and $VM,j=varj(SM,T)$⁠. Additionally, if we apply the property that the beta of the index portfolio is one, then $∑i=1Iwi,jβi,j=1$⁠, where $wi,j$ is the weight of security $i$ in the index at time $tj$⁠. This leads to Eq. (21).

Then, the covariance between security i and index M can be expressed as follows:

Finally, the realized third joint cumulant between the index M and an individual asset i in Eq. (8) can be expressed as follows:

So far, we have introduced methods to estimate realized higher-order (joint) cumulants using both option prices and stock prices. However, there is an issue with using the implied (co)moments discussed above. Prices are derived under the risk-neutral measure (Q-measure), whereas they evolve under the real measure (P-measure). Therefore, the implied (co)moments can differ from the true (co)moments under the P-measure. It also means that the practically obtained realized cumulants can differ from the true realized cumulants because they rely on Q-measure moments. In this regard, the next section investigates whether the realized cumulants obtained above can serve as good proxies for the true realized cumulants.

As discussed earlier, practically obtained realized cumulants can differ from the true realized cumulants. Therefore, this section investigates the validity of the realized cumulants obtained in practice through simulation. Details are as follows.

We assume the stochastic volatility model with contemporaneous jumps in volatility and prices (SVCJ), as proposed by Duffie et al. (2000), as follows:

Eqs. (24) and (26) represent the processes of the market index $SM,t$ and its instantaneous variance $vM,t$⁠, respectively. $ZM,t$ and $BM,t$ are Brownian motions with correlation $ρM$⁠; $NM,t$ follows a Poisson process with intensity $λM$⁠; $JM,t$ follows a normal distribution with mean $μM,S$ and standard deviation $σM,S$⁠; and $JM,v,t$ follows an exponential distribution with mean $μM,v$⁠. Eqs. (25) and (27) represent the stock price $Si,t$ and its instantaneous idiosyncratic variance $vi,t$⁠, respectively. $Zi,t$⁠, $Bi,t$⁠, $Ni,t$⁠, $Ji,t$⁠, and $Ji,v,t$ correspond to $ZM,t$⁠, $BM,t$⁠, $NM,t$⁠, $JM,t$⁠, and $JM,v,t$ of the market index. $βi,D$ and $βi,J$ represent the betas of diffusive risk and jump risk, respectively.

The parameters for the simulation are based on Broadie et al. (2007), Neuberger (2012), and Gourier (2016), as described in Table 1. Panel A shows parameters for a market index, and Panel B reports the parameters for idiosyncratic risk of two of the 29 individual stocks listed in Gourier (2016). Please refer to Gourier (2016) for the parameters of the remaining 27 stocks.

Using the P-measure parameters in Panel A and the first row (AXP) of Panel B in Table 1, we generate 10,000 paths of a pair consisting of an index and a security for 20 years on a daily basis. Then, we calculate and compare six estimators for each (joint) cumulant over the last month of each path. The estimators are as follows: Imp denotes the implied joint cumulant computed at the beginning of the last month. In practice, since the true model is unknown, one must extract this value from the prices of all Out-of-The-Money options with one month to maturity using the model-free approach of Bakshi et al. (2003). However, in the simulation we assume the SVCJ model, so we can obtain it directly by substituting the Q-measure parameters into Eq. (A1) in supplementary material. True denotes the true joint cumulant at the beginning of the last month. Although this value cannot be obtained in practice, it can be computed in the simulation by substituting the P-measure parameters into Eq. (A1). Conv-Real represents the conventional realized joint cumulants over the last month. They are obtained from the sub-period returns of the last month using the conventional comoments $∑j=1N(ΔSM,j)kM(ΔSi,j)ki$⁠. Hist represents the historical sample joint cumulants. They are estimated using monthly returns from the previous 24 months, measured at the beginning of the last month. Prac-Real represents the practical realized joint cumulants. They are estimated from daily returns and Imp, which is calculated for each day of the last month using Eq. (A1). Lastly, Real represents the realized (joint) cumulants. They are computed from daily returns and true moments for each day of the last month using Eq. (A1).

In Table 2, Panel A presents the mean and the standard deviation of the six estimators for the joint cumulants, and Panel B presents the corresponding values for coskewness and cokurtosis. In line with Eqs. (1)–(5), the average value of Real is close to True in Panel A ^[5]. On the other hand, each Conv-Real for the third joint cumulants is upwardly biased, whereas each Conv-Real for the fourth joint cumulant is downwardly biased. Additionally, every Conv-Real of coskewness is close to zero, and the Conv-Real cokurtosis is close to −3^[6]. These results show that Conv-Real fails to capture both contagion effects and volatility clustering effects, whereas True and Real do. Meanwhile, Prac-Real also deviates from True. However, the deviation of Prac-Real is smaller than that of Conv-Real or Imp. If we focus on the level of the estimators, Hist may be more accurate than Prac-Real. However, because Hist uses outdated data, it is not appropriate to represent the characteristics of realized return for a specific period when market conditions are time varying.

Table 3 shows the correlations among the (joint) cumulants reported in Table 2. Panel A compares the correlations between True and the other estimators of (co)skewness and (co)kurtosis. It shows that the correlation coefficients between True and Imp range from 0.941 to 0.993. Those between True and Hist are negative, ranging from −0.093 to −0.030. Those between True and Conv-Real range from −0.014 to 0.101. Those between True and Prac-Real range from 0.233 to 0.312. The correlations between True and Real range from 0.242 to 0.384. Therefore, Imp shows the highest correlation with True. This indicates that Imp and True move similarly, even though their levels differ, as shown in Table 2. This can be a natural result because True and Imp are ex ante estimators, whereas the others are ex post estimators obtained from realized price paths. Among the ex post estimators, Real is the most correlated with True. However, the correlations remain below 0.4. We conjecture that this is because Real reflects unexpected variations (especially jumps) in each period, whereas True does not. This may also explain the lower correlations between True and Real for kurtosis compared to skewness, since large values affect the higher-order moments more strongly.

According to the results in Panel A, Real represents the characteristics of realized returns for a period. However, Real is not empirically obtainable. To find an alternative to Real, Panel B of Table 3 shows the correlations between Real and the other ex post estimators. The correlation coefficients between Real and Prac-Real range from 0.797 to 0.940. By contrast, the correlation coefficient between Real and Hist (Conv-Real) is below 0.019 (0.195). Because Prac-Real is the most correlated with Real, we can use Prac-Real as a proxy for Real. Although they are estimated from different probability measures, the effect of using information under different measures does not appear to be severe.

To check the robustness of these results, we estimate (joint) cumulants using simulated prices with 29 different values of the idiosyncratic variation parameter from Panel B of Table 1. Panels C and D show the average values of 29 correlation coefficients among estimates. The results are consistent with Panels A and B. Specifically, in Panel C, Imp shows the highest correlation with True, while Prac-Real and Real exhibit higher correlation coefficients with True than Conv-Real. In Panel D, compared with Conv-Real and Hist, Prac-Real shows a higher average correlation with Real. These findings support the view that Prac-Real can be used as a proxy for Real and confirm that the effect of using implied moments instead of moments under the P-measure is limited.

In sum, True is an ex ante cumulant for a specific period, and it can be used for decision-making and planning. If market conditions are time-varying, Real is a useful ex post cumulant for a specific period, and it can be used for back-testing that period. Even if True and Real are not obtainable, we can use Imp and Prac-Real as their proxies, given the correlation observed earlier.

This section estimate realized (joint) cumulants using financial data and examine the characteristics of realized (joint) cumulants. First, we examine the realized cumulants of the S&P 500 index, estimated using daily index prices and S&P 500 European option data from January 1996 to August 2014. Option data, dividends, and the risk-free rate are obtained from OptionMetrics in Wharton Research Data Services (WRDS), and we construct a continuum of option prices using the methodology of Carr and Wu (2009) and Neuberger (2012).

In the analysis, we estimate various estimators of the cumulants of the S&P 500 index returns with three different time horizons: 30 days, 90 days, and 180 days. For example, for the cumulants of 30-day return: Imp is the implied cumulant extracted from option prices when the time to maturity of the options is 30 days; Conv-Real is estimated using daily index returns over the period; and Prac-Real is estimated using daily index returns together with option price data for the same period. Lastly, Sample represents the sample moments estimated using 30-day index returns over the whole period. As a result, for every 30-day window, we obtain the cumulants of Prac-Real, Conv-Real, and Imp, resulting in a total of 223 values for each cumulant over the whole period. By contrast, we have only one value for each cumulant in the case of Sample.

Panel A in Table 4 displays the averages of cumulants, skewness, and kurtosis for the 30-day return of the S&P 500. For the whole period, Sample skewness and kurtosis are −0.9 and 3.04, respectively. This indicates that the 30-day return distribution has historically been negatively skewed and leptokurtic, relative to the normal distribution. Similar results are observed for Prac-Real and Imp. The average implied skewness (Imp) and average Prac-Real skewness are −1.37 and −1.11, respectively, and the corresponding average kurtosis values are 5.94 and 4.33. Thus, the return distributions estimated by these two methods are also negatively skewed and leptokurtic. However, Conv-Real shows different results, as Conv-Real skewness and kurtosis are close to 0 and –3, respectively. This is consistent with the explanation in Footnote 6 and the results in Section 3. Panels B shows the averages of cumulants, skewness, and kurtosis of 90-day returns. Overall, all four estimators of the cumulants show similar patterns to those in Panel A. In other words, the return distribution is generally negatively skewed and leptokurtic, while the average Conv-Real skewness and kurtosis are close to 0 and –3, respectively. However, the return distribution becomes less negatively skewed and less leptokurtic as the time horizon increases, as reported in Campbell et al. (1997) and Munk (2013), consistent with the shape of the term structure of skewness and kurtosis in Neuberger and Payne (2021). Although not reported in the main text, this pattern persists when the analysis is conducted using 180-day returns. In addition, the standard deviations of kurtosis are larger than those of skewness for both Imp and Prac-Real. This is consistent with Schneider and Trojani (2019), because kurtosis is much more time-varying than skewness. This time-varying attribute is illustrated in Figure 1, which presents the variance, skewness, and kurtosis of 90-day returns as estimated by Imp and Prac-Real.

In addition, we examine whether Prac-Real can be explained by Imp and other realized moments as follows.

$Prac_Realt$ denotes the practical realized estimate (cumulant, joint cumulant, skewness, or kurtosis) at the end of month t. $Impt$ denotes the implied estimate at the beginning of month t. $Conv_Realt−1$ denotes the conventional realized estimate at the end of month t−1. $Histt$ denotes the sample estimate obtained at the beginning of month t using a sample of returns from the past 24 months.

The test above extends the analysis of realized variance in Jiang and Tian (2005) to the analysis of higher order estimator, and we could examine the relation between realized moments and implied moments and the predictability of realized moments. Kang and Lee (2016) and Kim (2016) show that realized higher-order moments are not predicted by other estimators. However, these results may not be accurate, since Conv-Real is a biased estimator and differs significantly from True, as shown in Section 3. We next examine the regression results of skewness and kurtosis for the S&P 500 index.

Table 5 reports the time-series regression results of Prac-Real skewness and kurtosis on Imp, lagged Prac-Real, lagged Conv-Real, and Hist. In univariate regressions, Imp and lagged Prac-Real have predictive power for the Prac-Real skewness, unlike the regression results for Conv-Real in previous studies. In multivariate regressions, Imp and lagged Prac-Real still have predictive power for the skewness of future returns, consistent with Neuberger (2012). For kurtosis, in univariate regressions, Prac-Real is explained by Imp and lagged Prac-Real at the 5% significance level, and by lagged Conv-Real and Hist at the 10% significance level. In multivariate regressions, Prac-Real is predicted by Imp and lagged Prac-Real, consistent with the results for skewness. In sum, implied and lagged realized skewness (kurtosis) have predictive power for realized skewness (kurtosis), unlike the findings of previous studies.