Approximate density from moments and quantiles, then sample from it

Question

Situation

I need to send R code to a third party to run estimations for me (I will not be able to work with the data directly). I want to simulate data to test some of the estimators before sending the code to them.

The data provider has given me the following summary statistics for all variables: The first four moments (mean, variance, skewness, kurtosis), the four largest/smallest values (min/max), and the 1%, 5%, 10%, 25%, 50%, 75%, 90%, 95%, 99% percentiles.

I believe the output was obtained in Stata, I am pasting an example below. The example is not from the actual data (actual data has close to 6 million observations).

Question

What would be the best way to simulate data? Initially I was just going to pick a distribution and sample from it (e.g. binomial/multinomial/(log-)normal/exponential/truncated normal), but with the information provided I presumed it is possible to do better, at least for variables which aren't binomial.

Input example

      Percentiles      Smallest
 1%          163             99
 5%          216            111
10%          248            113       Obs               3,170
25%          322            114       Sum of Wgt.       3,170

50%          494                      Mean           1262.359
                        Largest       Std. Dev.      3093.165
75%          984          41584
90%         2450          54413       Variance        9567670
95%         5181          58477       Skewness       10.59025
99%        10826          59349       Kurtosis       157.7004

What I am currently doing

I have tried fitting the parameters of a specific distribution I think might be appropriate using functions in the library(rriskDistributions) (e.g. get.lnorm.par()). This works well sometimes but often it does not.

Currently I am fitting the CDF using splines, obtain the PDF using the spline functions derivative, and then sample from it.

Neither of these approaches work very well in general. I'm hoping for an approach that is generic and delivers a good approximation without me manually eyeballing the distribution and investigating the fit. I understand that this may be asking a lot given the limited data at my disposal.

## function for spline interpolation
splsample <- function(p, v,
                      size = 1000000,
                      vmin = min(v), vmax = max(v),
                      gridsize = min(3*(vmax-vmin), 1000),
                      step = NULL, plot = FALSE, ...) {
  s <- splinefun(v, p, ...)
  if(is.null(step)) {
    grid <- seq(from = vmin, to = vmax, length.out = gridsize)
  } else {
    grid <- seq(from = vmin, to = vmax, by = step)
  }
  pr <- s(grid, deriv = 1)
  pr[pr < 0] <- 0
  if (plot == TRUE) {
    plot(grid, pr)
  }
  bs <- sample(grid, p = pr, size = size, replace = TRUE)
  return(bs)
}

## input
percentiles <- c(0.01, 0.05, 0.10, 0.25, 0.5, 0.75, 0.9, 0.95, 0.99)
values <- c(163, 216, 248, 322, 494, 984, 2450, 5181, 10826)

## spline approximation of pdf
x <- splsample(percentiles, values, plot = TRUE)
summary(x)
mean(x)
var(x)

## alternative: fitting a truncated normal
library("rriskDistributions")
library("msm")
dpar <- get.tnorm.par(p = percentiles, q = values)
x <- rtnorm(10000, mean = dpar["mean"], sd = dpar["sd"],
            lower = dpar["lower"], upper = dpar["upper"])
x[x < 0] <- 0
summary(x)
mean(x)
var(x)

Answer 1

To quickly simulate based on moments, try rpearson() from library(PearsonDS).

library(PearsonDS)
target.moms <- c(1262.39, 9567670, 10.59025, 157.7004)
y <- rpearson(n=1000000, moments=target.moms)

rpearson() works well for matching the moments. However, the splines approach that you're already using will be better at recovering the percentiles. See below for an example.

#Evaluating the results
library(moments)
eval <- function(data) {
  result.list <- list(mean=mean(data),
                   var=var(data),
                   skew=skewness(data),
                   kurt=kurtosis(data),
                   quantile(data, c(.01,.05,.10,.25,.50,.75,.90,.95,.99) ) )
  round(unlist(result.list), 2)                
}
x <- splsample(percentiles, values, plot = TRUE)
eval(x) #splines
eval(y) #rpearson()