Unbiased weighted variance and skewness estimators for overlapping returns

This article develops unbiased weighted variance and skewness estimators for overlapping return distributions. These estimators extend the variance estimation methods constructed in Bod et. al. (Applied Financial Economics 12:155-158, 2002) and Lo and MacKinlay (Review of Financial Studies 1:41-66, 1988). In addition, they may be used in overlapping return variance or skewness ratio tests as in Charles and Darné (Journal of Economic Surveys 3:503-527, 2009) and Wong (Cardiff Economics Working Papers, 2016). An example using synthetic overlapping returns from a model fit to data from the SPY S&P 500 exchange traded fund is given in order to demonstrate under which circumstances the unbiased correction becomes significant in skewness estimation. Finally, we compare the effect of the HAC weighting schemes of Andrews (Econometrica 53:817-858, 1991) as a function of sample size and overlapping return window length.


Introduction
Overlapping returns are used in many contexts in the finance and econometrics literature. Applications include variance ratio tests, regression parameter error estimation, and alternative resampling methods. Standard statistical inference and estimation techniques applied to overlapping return financial time series are typically biased. In addition, for such series, recent data is regularly viewed as more relevant than past information, which has resulted in the creation of weighted generalizations of estimation methodologies. This motivates the development of unbiased analogues of such estimators which we explore in the cases of the variance and skewness statistics. Our central aim is to construct unbiased weighted variance and skewness estimators for overlapping return distributions.
Several estimation procedures and hypothesis testing frameworks have been improved through the utilization of overlapping returns. In financial overlapping return applications, Lo and MacKinlay (1988) and Hansen and Hodrick (1980) demonstrate how overlapping returns may be used to increase the efficiency of statistics used in *Correspondence: smt@njit.edu Martin Tuchman School of Management at the New Jersey Institute of Technology, Newark, NJ, USA variance ratio tests. Dunis and Keller (1995) developed a panel regression method based on overlapping returns, and Müller (1993) concludes utilizing overlapping returns in most applications will result in an overall increase in estimation precision of statistics that are a function of the overlapping returns when compared with their analogues for simple returns. In Jackwerth (2000), the author discovered that the overlapping return distribution for the S&P 500 is left-skewed and examined differences between risk neutral and realized distributions between overlapping and non-overlapping returns of the S&P 500 index and observed how associated risk aversion functions changed dramatically around the 1987 stock market crash. Wong (2016) develops skewness and kurtosis ratio tests for overlapping returns. The new weighted unbiased skewness estimator constructed below may be used as an input into any of these applications.
The idea of assigning greater weight to recent data and less weight to past data has been discussed in a number of econometric and financial studies. Past economic data may have little impact or be entirely irrelevant for present projections. In addition, by placing additional weight on recent data, associated estimation procedures tend to react more strongly to structural changes in the underlying assumption about the distribution the sample is drawn from than their uniformly weighted counterparts. For example, Tsokos (2010) shows that under nonstationary economic realization, weighted moving average models perform significantly better than the classical ARIMA model in forecasting stock prices. Andrews (1991) develops weighting schemes used in the estimation of covariance matrices assuming the underlying time series exhibits nontrivial autocorrelation and heteroskedasticity which we will utilize below. Weighted estimators are routinely used in practice as well. In particular, in Longerstaey and Spencer (1996), it is demonstrated that exponentially weighted moving average estimators incorporate external shocks more readily than equally weighted moving averages, thus providing a more realistic measure of current volatility.
Volatility and skewness estimation of financial return distributions has been the subject of a number of articles. Early examples include using maximum likelihood estimation to fit a model distribution to observed data and computing the associated model statistics in Fama (1965) and Mandelbrot (1963). More recent work has focused on the estimation of stochastic volatility models in Broto (2004). Time series techniques have also been widely applied to this task, c.f. (Tsay 2010). Measuring the asymmetry of financial return distributions has also been the central theme of many references. Grigoletto and Lisi (2006) and Wen and Yang (2009) find persistent non-trivial skewness is present in the simple daily return distributions of nearly every major international equity index. Xu (2007) shows that equity return distribution skewness is positively correlated with simultaneous returns and negatively correlated with lagged returns.
When working with overlapping returns, especially when encountering small sample sizes, bias effects from standard estimators, such as the sample variance, become important. In Lo and MacKinlay (1988), the authors provide a consistent but biased overlapping return variance estimator that has been used in several subsequent references, including Liu and He (1991), Fong et al. (1997), and Amélie and Olivier (2009). This estimator was improved in Bod et al. (2002) where the authors constructed an unbiased variance estimator for unweighted overlapping returns. Kluitman and Franses (2002) extended this work to develop an estimator that includes the case where the returns have nontrivial autocorrelation. Our main contribution is to extend these results by developing weighted unbiased variance and skewness estimators for overlapping return time series.
This article is organized as follows. We first fix notation and then derive an unbiased weighted estimator for the variance of a time series of overlapping returns. We give reduced expressions for this estimator in the cases of uniform and exponential weights. Next, we construct a similar weighted unbiased estimator for the skewness of an overlapping return distribution. We then demonstrate the difference between a normalized version of the skewness estimator and the standard normalized sample skewness in a simulation which models the overlapping return distribution of the S&P 500 index, and then summarize our results. We finally compare the estimation of the weighted volatility and skewness of the overlapping return distribution of the S&P 500 index for various weighting schemes, sample sizes, and overlapping lengths and conclude with potential additional questions to explore.

Methodology
We begin by establishing notation. Given integers n, q > 0 with q < n, let p t > 0 for t = 0, . . . , n + q − 1 denote an asset price time series by p t and let r t = p t /p t−1 − 1 be the associated simple returns where t = 1, . . . , n + q − 1. Following Bod et al. (2002) and Lo and MacKinlay (1988), we assume that r t have zero mean, E[r t ] = 0, covariance E[r t r s ] = 0 for any t > s, and equal finite variance Var(r t ) = σ 2 < ∞. We note that in Lo and MacKinlay (1988) and subsequent references, the authors show that these assumptions, referred to as the random walk version of the martingale hypothesis, do not hold for a variety of financial time series. This is achieved by assuming a null hypothesis that they hold and then demonstrating how variance ratios of overlapping returns may be used to reject this assertion. In this spirit, we proceed by defining the q-period overlapping returns y t associated with r s by ( 1 ) We construct weighted unbiased estimators of the variance and skewness of y t and pair a weight w t with each y t such that w t > 0 and n t=1 w t = 1. Let W ts = s k=t w k be the sum of the t-th through s-th weight, and note W 1n = 1.
We first derive an unbiased weighted overlapping return variance estimatorσ 2 y . We seek an estimator of the form where we find C 1 such thatσ 2 y is an unbiased estimator of qσ 2 . Note that since r i are independent, the true variance of each y t is Var(y t ) = qσ 2 so thatσ 2 is an unbiased estimator of the variance of the random variable from which y t are sampled from if This may be viewed as a constraint that defines the constant C 1 (n, q, w) which we determine by computing and noting that Combining the above, we find where we note that E(y t ) = E (ȳ w ) = 0, and the last equality follows from E y 2 t = Var(y t ) = qσ 2 as well as In order to compute the variance term in the previous equation, we decompose the weighted average of y t as One may arrive at this decomposition by viewing the individual return terms of the sumȳ w = t w t y t as a table with values w t r s whose first row has elements, w 1 r 1 , w 1 r 2 , . . . , w 1 r q and final row is given by w n r n , w n r n+1 , . . . , w n r n+q−1 . Note thatȳ w is equivalent to the sum of all values in this table. The first term in this decomposition corresponds to grouping all elements of this table above the diagonal whose edge is formed by the w 1 r q and w q r q entries and factoring our common returns multiplied into varying weights. The final term can be arrived at by aggregating all terms below the diagonal formed by the w n−q+2 r n+1 and w n r n+1 entries. The middle sum is computed by combining the remaining terms in the table.
Since each individual sum is composed of returns that are independent of all the returns in the other two sums, we find that Solving for the unbiased constant from Eqs. (7) and (9), we arrive at We note that in the case of uniform weights w t = 1/n, this result reduces to that of Bod et. al. (2002), where We now derive an unbiased skewness estimator in a similar manner. Assume that E r 3 t = γ , so that E y 3 t = qγ , and consider an estimator of the skewness of the overlapping return distribution of the form We would like to construct an unbiased estimator of qγ which requires To determine C 2 , we compute Decomposing y t into independent r t sums allows one to calculate the first term with For the second term, we note that Computing the expectation, we find, Note that only terms of the form r 3 t are non-trivial in expectation. Splitting into two cases, we have Combining Eqs. (16) and (20), we find that For the final two terms, note E y t ȳ w 2 = s,k w s w k E(y t y s y k ), E ȳ w 3 = t,s,k w t w s w k E(y t y s y k ), so we may simplify n t=1 w t 3y t ȳ w 2 − ȳ w 3 = 2E ȳ w 3 .
Using the decomposition in Eq. (14), by independence of the terms Combining results from Eqs. (15), (21), and (24), we arrive at an expression for C 2 given by C 2 (n, q, w) = 1 − 3 q n t,s=1 |t−s|<q w t w s (q − |t − s|) In the case of uniform weights w t = 1/n, the terms in this expression simplify to n t,s=1 |t−s|<q which yields a simple form for C 2 given by We finally note that it is possible to derive a closed form expression for C 2 in the case of exponential weights that were previously considered for the variance estimator; however, the expression is quite lengthy and we omit it here but it is available upon request. Finally, we turn to a simulation in order to understand for what parameter pairs (n, q) the effects of the unbiased skewness estimator are most significant.

Empirical studies and results
We now develop a simulation to compare the relative error of the uniformly weighted unbiased skewness estimator γ y and the standard unbiased sample skewness estimator which may be found in Zwillinger and Kokoska (2000). We first construct a dataset of end of day simple returns calculated from closing prices for the SPY exchange traded fund from January 1, 2012, to December 31, 2016. This was achieved using Bloomberg's Python API and an associated wrapper package named tia. We downloaded historical end of day closing prices identified with Bloomberg's PX_LAST field that are both split and dividend adjusted. This time series was fully populated with data, and hence, there was no need to fill in missing values.
Next, we fit a model distribution to this data in order to establish a framework to test the weighted unbiased skewness and kurtosis estimators in a setting where the true values of these statistics are known which closely approximates actual market data. To this end, we let X denote a skew normal distribution whose probability density function is given by where here φ and are the probability and cumulative distribution functions of a standard normal random variable and b > 0. The skew normal distribution has mean μ and variance σ 2 defined by μ = a + bd 2 π , The normalized skewness γ of this distribution is given by We fit this distribution to the SPY daily return data using maximum likelihood estimation. Let the likelihood function of this model be denoted by L. We find the model parameters by directly maximizing the log-likelihood function using a BFGS optimizer developed in Fletcher (1987) and determine that the best fit parameters are given byâ = 0.00640,b = 0.01006, andĉ = −1.1029, which have associated mean, variance, and normalized skewness given by 4.528 × 10 −4 , 6.5789 × 10 −5 , and − 0.1689, respectively. Returns will be drawn from this distribution in the Monte Carlo simulation considered below. Specifically, this simulation consists of sampling n values from the MLE distribution, computing the q-period overlapping returns of these time series, then calculating the normalized unbiased sample skewnessγ s (see Zwillinger and Kokoska (2000)) and the normalized unbiased skewnessγ y with uniform weights given byγ y =γ y /σ 3 y . We then compute the relative error |γ s −γ y |/γ y , between the different estimators for each simulation, repeat the above process 10,000 times, and display the mean percentage errors in Table 1 for distinct (n, q) pairs.
We omit cases where the number of overlapping returns n − q ≤ q, and first note that as the sample size increases, the error between the two estimators decreases for any fixed q value. However, when q/n is relatively large, say greater than 5%, then there are significant differences between the two estimators.
Next, we explore several weighting schemes described in Andrews (1991) which are widely used for covariance matrix estimation in the presence of heteroskedasticity and autocorrelation. Specifically, we consider weights constructed by Bartlett (Oppenheim et al. 1999), Parzen (White 1980), Tukey-Hamming (Blackman and Tukey 1958), and the Quadratic Spectral weights of Priestley (1962) and Epanechnikov (1969). These weights are defined in terms of a kernel function k(·) and are given by w t = k(bt/T) where T is a bandwidth parameter and b is a scaling constant in Zeileis (2004) and Zwillinger (2000). There are many references that study the problem of optimal bandwidth selection c.f. (Lazarus et al. 2017;Newey and West 1994;Stock and Watson 2011;Wooldridge 2006); however, we are interested in constructing reasonable weighting schemes to place importance on more recent over prior data. We found that setting the bandwidth to the sample size and b = 1.2 achieves this aim.
In Fig. 1, we plot the kernel functions associated with these weighting schemes given in Andrews (1991) which are defined by k QS (x) = 25 12π 2 x 2 sin(6πx/5) 6πx/5 − cos(6πx/5) . (37) Note that when using weighted estimators, one effectively reduces the original sample size. For example, in the extreme case of binary zero or one valued weights, only the weights with value one contribute to the estimator which reduces the sample size to the percentage of one valued weights. In this example, we can find the percentage sample size reduction by approximating the area under each weight curve. In reference to the uniform weights which we take to have normalized area of 1, the HAC weights effectively reduce the sample size by PR: 31%, BT,TK: 41%, and QS: 54% so that uniform weights have approximately two to three times the sample size of these weighting schemes. Next, we examine how estimation of the unbiased weighted standard deviation and skewness estimators varies as a function of the overlapping return period q, sample size n, and weighting scheme using the SPY dataset previously described. We consider overlapping return periods of 5, 21, and 63 sample points which correspond to weekly, monthly, and quarterly aggregation windows for our example daily return data. Next, we truncate the sample size to 256, 512, and 1024 data points which roughly correspond to 1, 2, and 5-year time periods, using a trailing truncation window on the SPY returns.
Estimation results are presented in Table 2 where overlapping return standard deviations are displayed as percentages. We first note that estimates for uniform weights tend to be outliers as they include two to three times the effective sample size as the other weighting schemes. Next, note that overall standard deviation values are greatest in the n = 512 case, slightly lower for n = 256, and considerably lower for n = 1024. This is due to the historical volatility of the S&P 500 being relatively high during the mid 2015 to early 2016 time period and lower in prior years over the five year window being considered. One may also note that the overlapping return distribution becomes increasingly negatively skewed as a result.
Finally, we compare how skewness estimation varies over time for different weighting schemes. In particular, we consider a 4-year window from 1/1/2012 to 1/1/2016 and estimate the unbiased skewness using a 252-day lookback period for each of the HAC estimators and a 126-day lookback period for the uniformly weighted estimator which is selected to ensure that the sample sizes are on par with one another. This procedure is carried out on a rolling basis, and results are plotted in Fig. 2.
We note that the general forms of the time series in Fig. 2 tend to be similar for the majority of dates displayed. The HAC weighted estimators are more reactive than the uniform estimator and do not exhibit single day jumps that are as large in magnitude as the HAC uniform estimator.

Discussion and conclusions
In summary, we have derived closed form expressions for weighted unbiased variance and skewness estimators. We also developed simplified expressions for these estimators in the case of exponential weights for the variance estimator and uniform weights for both estimators. The differences between the standard unbiased sample skewness and new normalized unbiased skewness estimators were demonstrated to be significant in the case of skewness estimation for SPY end of day return data for HAC weighting schemes.
We note that as in Bod et al. (2002) and Lo and MacKinlay (1988), we assume returns satisfy the random walk version of the martingale hypothesis which generally does not hold for financial time series. An interesting future application of the skewness estimator would be to develop a hypothesis test for this assumption which may compliment the results in Lo and MacKinlay (1988). For additional future work, it would be of interest to consider as in Kluitman and Franses (2002) analogues of the weighted unbiased variance and skewness estimators under the assumption that the return process satisfies an AR(1), MA(1) or more general time series model. This will require repeating t he above derivations and retaining terms of the form E [r t r s ] and E [r t r s r k ] for t = s = k which are no longer trivial but depend on the underlying model one assumes for the return process. Then, one could fit such models to market data and compare values of the two estimators. We anticipate the estimator will strongly depend on the sign of the AR(1) lag parameter and the Fig. 2 Rolling n = 252 day overlapping return unbiased skewness for HAC weighting schemes and a n = 126 lookback period for uniform weights with a q = 21 overlapping period for SPY daily returns white noise parameter of the MA(1) model as shown in Kluitman and Franses (2002) but leave this for a future study.