Gilgamath

I recently pushed version 0.2.0 of my fromo package to CRAN. This package implements (relatively) fast, numerically robust computation of moments via Rcpp.

The big changes in this release are:

• Support for weighted moment estimation.
• Computation of running moments over windows defined by time (or some other increasing index), rather than vector index.
• Some modest improvements in speed for the 'dangerous' use cases (no checking for NA, no weights, etc.)

The time-based running moments are supported via the t_running_* operations, and we support means, standard deviation, skew, kurtosis, centered and standardized moments and cumulants, z-score, Sharpe, and t-stat. The idea is that your observations are associated with some increasing index, which you can think of as the observation time, and you wish to compute moments over a fixed time window. To bloat the API, the times from which you 'look back' can optionally be something other than the time indices of the input, so the input and output size can be different.

Some example uses might be:

• Compute the volatility of an asset's returns over the previous 6 months, on every trade day.
• Compute the total monthly sales of a company at month ends.

Because the API also allows you to use weights as implicit time deltas, you can also do weird and unadvisable things like compute the Sharpe of an asset over the last 1 million shares traded.

Speed improvements come from my random walk through c++ design idioms. I also implemented a 'swap' procedure for the running standard deviation which incorporates a Welford's method addition and removal into a single step. I do not believe that Welford's method is the fastest algorithm for a summarizing moment computation: probably a two pass solution to compute the mean first, then the centered moments is faster. However, for the …

Another confidence limit on the Signal Noise ratio of the Markowitz portfolio.

In a previous blog post, I looked at asymptotic confidence intervals for the Signal to Noise ratio of the (sample) Markowitz portfolio, finding them to be deficient. (Perhaps they are useful if one has hundreds of thousands of days of data, but are otherwise awful.) Those confidence intervals came from revision four of my paper on the Asymptotic distribution of the Markowitz Portfolio. In that same update, I also describe, albeit in an obfuscated form, the asymptotic distribution of the sample Markowitz portfolio for elliptical returns. Here I check that finding empirically.

Suppose you observe a \(p\) vector of returns drawn from an elliptical distribution with mean \(\mu\), covariance \(\Sigma\) and 'kurtosis factor', \(\kappa\). Three times the kurtosis factor is the kurtosis of marginals under this assumed model. It takes value \(1\) for a multivariate normal. This model of returns is slightly more realistic than multivariate normal, but does not allow for skewness of asset returns, which seems unrealistic.

Nonetheless, let \(\hat{\nu}\) be the Markowitz portfolio built on a sample of \(n\) days of independent returns:

where \(\hat{\mu}, \hat{\Sigma}\) are the regular 'vanilla' estimates of mean and covariance. The vector \(\hat{\nu}\) is, in a sense, over-corrected, and we need to cancel out a square root of \(\Sigma\) (the population value). So we will consider the distribution of \(Q \Sigma^{\top/2} \hat{\nu}\), where \(\Sigma^{\top/2}\) is the upper triangular Cholesky factor of \(\Sigma\), and where \(Q\) is an orthogonal matrix (\(Q Q^{\top} = I\)), and where \(Q\) rotates \(\Sigma^{-1/2}\mu\) onto \(e_1\), the first basis vector:

where \(\zeta\) is the Signal to Noise ratio of the population Markowitz portfolio: \(\zeta = \sqrt{\mu^{\top}\Sigma^{-1}\mu} = \left\Vert …

I wanted a drop-in replacement for geom_errorbar in ggplot2 that would plot a density cloud of uncertainty. The idea is that typically (well, where I work), the ymin and ymax of an errorbar are plotted at plus and minus one standard deviation. A 'cloud' where the alpha is proportional to a normal density with the same standard deviations could show the same information on a plot with a little less clutter. I found out how to do this with a very ugly function, but wanted to do it the 'right' way by spawning my own geom. So the geom_cloud.

After looking at a bunch of other ggplot2 extensions, some amount of tinkering and hair-pulling, and we have the following code. The first part just computes standard deviations which are equally spaced in normal density. This is then used to create a list of geom_ribbon with equal alpha, but the right size. A little trickery is used to get the scales right. There are three parameters: the steps, which control how many ribbons are drawn. The default value is a little conservative. A larger value, like 15, gives very smooth clouds. The se_mult is the number of standard deviations that the ymax and ymin are plotted at, defaulting to 1 here. If you plot your errorbars at 2 standard errors, change this to 2. The max_alpha is the alpha at the maximal density, i.e. around y.

# get points equally spaced in density 
equal_ses <- function(steps) {
    xend <- c(0,4)
    endpnts <- dnorm(xend)
# perhaps use ppoints instead?
    deql <- seq(from=endpnts[1],to=endpnts[2],length.out=steps+1)
    davg <- (deql[-1] + deql[-length(deql)])/2
# invert
    xeql <- unlist(lapply(davg,function(d) {
                     uniroot(f=function(x) { dnorm(x) - d },interval=xend)$root
    }))
    xeql
}

library(ggplot2)
library(grid)

geom_cloud <- function(mapping …