Structural analysis of first-price auctions with qdfauction: the USFS timber data

This vignette walks through the full pipeline of qdfauction on real data: the sealed-bid subsample of the U.S. Forest Service timber auctions assembled by Lu and Perrigne (2008). We recover bidders’ private values from their bids, estimate the private-value density, and estimate the seller’s expected revenue as a function of the reserve price.

The data contain 107 auctions with two bidders and 108 with three. We analyse the two groups separately so that the number of bidders is known and fixed within each group. The bids are homogenised: within each group the log bids were regressed on the log appraised per-unit value and the log timber volume of the lot, and rescaled to the average lot, so that the bids of a group can be treated as draws from a common distribution. See ?timber.

1. The quantile density of bids

In a symmetric first-price auction with independent private values the inverse bid function is (Guerre, Perrigne and Vuong 2000)

\[v = b + \frac{1}{n-1}\frac{G(b)}{g(b)} = b + \frac{u\,q(u)}{n-1},\qquad u = G(b),\]

where \(q = 1/g(Q(u))\) is the quantile density of bids. The second form (Zhang 2022; Doosti, Dewan and Talebian 2025) avoids dividing by a density estimate that is small and unstable near the boundaries. qdf() estimates \(q\) with ten estimators; the indirect Poisson estimator with biased cross-validation performed best in the 2025 simulations.

b3  <- as.numeric(timber$n3)
# bandwidth = "bcv" selects h by biased cross-validation (about 20 s here);
# we fix h at the value it selects to keep the vignette fast
q_ip <- qdf(b3, method = "indirect_poisson", bandwidth = 0.02139)
q_b  <- qdf(b3, method = "bernstein",        bandwidth = "bcv")
q_w  <- qdf(b3, method = "wavelet_block",    bandwidth = 0.15)
q_ip
#> Quantile density estimate
#>   method    : indirect_poisson 
#>   n         : 324 
#>   scale     : 55.57 (estimated on x / scale)
#>   h (1/m)   : 0.02139  (m = 46.75) 
#>   selector  : fixed bandwidth
#>   q(u) on 324 grid points; range     61.15 to 144079.65
plot(q_ip, ylim = c(0, 400), main = "Quantile density of bids, N = 3")
lines(q_b$u, q_b$q, col = 2); lines(q_w$u, q_w$q, col = 4)
legend("topleft", c("indirect Poisson", "Bernstein", "block wavelet"),
       col = c(1, 2, 4), lty = 1, bty = "n")

The quantile density rises steeply toward \(u = 1\): the bid density is close to zero at the top of the support, which is exactly where the conventional ratio \(G/g\) is unstable.

2. Recovering private values

fpa_values() applies the inversion. We compare the indirect Poisson inversion with the conventional GPV estimator and with the boundary-corrected kernel of Hickman and Hubbard (2015), and impose monotonicity of the inverse bid function with the natural-spline fix of the 2025 paper.

fv_ip  <- fpa_values(timber$n3, method = "indirect_poisson", bandwidth = q_ip$h,
                     monotone = "spline")           # reuse the BCV bandwidth from above
fv_gpv <- fpa_values(timber$n3, method = "gpv")
fv_ip
#> Private values recovered from first-price auction bids
#>   bids      : 324 pooled, 3 bidders per auction
#>   method    : indirect_poisson (h = 0.02139) 
#>   utility   : risk neutral 
#>   copula    : independence 
#>   monotone  : spline 
#>   markup    : median 16.81 ( 25.4 % of bid )
plot(fv_ip, main = "Inverse bid function, N = 3")
points(fv_gpv$bids, fv_gpv$values, col = 2, pch = 3, cex = 0.6)
legend("topleft", c("indirect Poisson (raw)", "monotone fit", "GPV"),
       col = c("grey50", 4, 2), pch = c(20, NA, 3), lty = c(NA, 1, NA), bty = "n")

The implied markups \(\hat v - b\) fall with competition, as theory predicts:

fv2 <- fpa_values(timber$n2, method = "indirect_poisson", bandwidth = 0.02108,   # BCV choice
                  monotone = "spline")
rbind(`N = 2` = quantile(fv2$values - fv2$bids, c(.25, .5, .75)),
      `N = 3` = quantile(fv_ip$values - fv_ip$bids, c(.25, .5, .75)))
#>             25%      50%      75%
#> N = 2 15.940769 34.15681 93.37416
#> N = 3  9.397149 16.81135 47.16329

Lu and Perrigne (2008) find evidence of risk aversion in these auctions. Under CRRA utility \(U(t) = t^{1-\eta}\) the inversion becomes \(v = b + (1-\eta)\,u q(u)/(n-1)\), so risk-neutral pseudo-values are upper bounds on the true ones. The eta argument handles this:

fv_ra <- fpa_values(timber$n3, method = "indirect_poisson", bandwidth = q_ip$h,
                    eta = 0.3, monotone = "spline")
median(fv_ra$values - fv_ra$bids) / median(fv_ip$values - fv_ip$bids)   # = 1 - eta
#> [1] 0.7

3. The private-value density

fpa_density() estimates the density of private values in two steps. The first step recovers pseudo-values (Hickman–Hubbard boundary-corrected GPV by default, or any fpa_values() fit); the second step estimates their density with the wavelet estimators of Doosti, Dewan and Sbaï (2026) or one of the kernel competitors. Because the bids are in dollars, they are mapped to \([0,1]\) by a log map and the density is returned in dollars through the Jacobian.

fd_w  <- fpa_density(timber$n3, method = "wavelet_ll", j0 = 4, h = 0.15)
fd_hh <- fpa_density(timber$n3, method = "kde_hh")
fd_ms <- fpa_density(timber$n3, method = "marmer_shneyerov")
fd_q  <- fpa_density(timber$n3, method = "wavelet_ll", first_stage = fv_ip, j0 = 4, h = 0.15)
fd_w
#> Private-value density, first-price auction
#>   bids        : 324 pooled, 3 bidders per auction
#>   first stage : gpv_hh 
#>   method      : wavelet_ll (j0 = 4, h = 0.15) 
#>   scale map   : log 
#>   pseudo-values:    7.639 8645.868 ; max v / max b = 18.2
plot(fd_w, xlim = c(0, 400), ylim = c(0, 0.02), main = "Private-value density, N = 3")
#> Warning in graphics::rug(x$values, col = "grey60"): some values will be clipped
lines(fd_hh$x, fd_hh$f, col = 2, lty = 2)
lines(fd_ms$x, fd_ms$f, col = 3, lty = 3)
lines(fd_q$x,  fd_q$f,  col = 4, lty = 4)
legend("topright", c("wavelet LL (GPV-HH first stage)", "GPV-HH kernel",
                     "Marmer-Shneyerov", "wavelet LL (QDF first stage)"),
       col = 1:4, lty = 1:4, bty = "n")

The tuning parameters of the wavelet estimator can be chosen by cross-validated held-out log-likelihood (j0 = "cv"), and the estimators can be compared on the same footing with fpa_density_cv():

fd_cv <- fpa_density(timber$n3, method = "wavelet_ll", j0 = "cv")
fpa_density_cv(timber$n3, methods = c("wavelet_ll", "kde_hh", "marmer_shneyerov"),
               folds = 5, j0 = fd_cv$j0, h = fd_cv$h)

Confidence bands

For the GPV estimator, fpa_density_ci() adds analytic standard errors, pointwise confidence intervals and a bootstrap uniform confidence band (Ma, Marmer and Shneyerov 2019). The band is valid over an interior range of values; by default the grid spans the 10th to 90th percentiles of the trimmed pseudo-values.

ci3 <- fpa_density_ci(timber$n3, B = 199)
ci3
#> GPV estimate of the private-value density with confidence bands
#>   (Ma, Marmer and Shneyerov 2019)
#>   bids        : 324 pooled, 3 bidders per auction
#>   bandwidths  : h_g = 29.76  h_f = 49.96 
#>   grid        : 49.82 to 227.3 ( 201 points )
#>   level       : 0.95 ; pointwise intervals: percentile 
#>   uniform band: critical value 3.252 from 199 bootstrap samples
#>                 max width 0.005014
plot(ci3, main = "GPV density with 95% uniform band, N = 3")

4. Expected revenue and the optimal reserve price

Given recovered values, fpa_revenue() estimates the seller’s expected revenue as a function of the reserve price (Zincenko 2024; Doosti 2026b), and fpa_optimal_reserve() maximises it. Here we take the seller’s own valuation to be v_s = 20 dollars per unit.

r  <- seq(20, 120, by = 5)
R3 <- fpa_revenue(fv_ip, r, v_s = 20)
opt <- fpa_optimal_reserve(fv_ip, v_s = 20)
plot(r, R3, type = "l", lwd = 2, xlab = "reserve price ($/unit)",
     ylab = "expected revenue", main = "Expected revenue, N = 3")
abline(v = opt$r_star, lty = 2)

opt[c("r_star", "revenue")]
#> $r_star
#> [1] 33.65819
#> 
#> $revenue
#> [1] 93.32423

Bootstrap standard errors, pointwise intervals and a uniform band follow from fpa_revenue_boot(); the Bernstein estimator is used here because it is fast enough for many bootstrap replications. The bandwidth is selected once, by BCV on the full sample, and held fixed in every bootstrap replication. If the criterion is minimised at an end of m_range, fpa_revenue_boot() warns and records it: such an m was not really selected by the data.

rb <- fpa_revenue_boot(timber$n3, r = seq(30, 100, by = 10), B = 30, v_s = 20,
                       method = "bernstein", bandwidth = "bcv")
rb$bandwidth[c("m", "at_boundary", "boundary")]
#> $m
#> [1] 40.75712
#> 
#> $at_boundary
#> [1] FALSE
#> 
#> $boundary
#> [1] NA
rb
#> Expected revenue with bootstrap inference (B = 30 , resampling auctions )
#>   bandwidth    : m = 40.76 selected by bcv on the full sample, fixed in all replications
#>   pointwise CI and uniform band at level 0.95
#>       r estimate     se ci_lower ci_upper band_lower band_upper
#> 30   30  97.1882 7.7523  85.2409 109.1355    82.6128   111.7636
#> 40   40  97.0085 7.7196  85.3131 108.7040    82.4944   111.5226
#> 50   50  96.3438 7.9530  84.8110 107.8766    81.3910   111.2966
#> 60   60  96.2322 7.8205  84.9189 107.5456    81.5286   110.9359
#> 70   70  91.7032 8.1009  79.4719 103.9346    76.4724   106.9341
#> 80   80  91.2244 8.5870  76.2672 106.1815    75.0794   107.3693
#> 90   90  87.1835 9.2213  70.9909 103.3761    69.8461   104.5210
#> 100 100  85.8091 9.1938  68.6358 102.9824    68.5233   103.0948
plot(rb, main = "Expected revenue with 95% pointwise CI and uniform band")

5. Checking the method by simulation

fpa_simulate() draws values from any marginal (and, optionally, an Archimedean copula for affiliated values or a CRRA risk-aversion parameter) and computes exact equilibrium bids, and fpa_montecarlo() scores a list of estimators over replications. The following reproduces one row of Table 2 of Doosti, Dewan and Talebian (2025) in a few seconds:

pdf <- function(x) dgamma(x, 10, 2); cdf <- function(x) pgamma(x, 10, 2)
mc <- fpa_montecarlo(R = 3, L = 20, n = 5, Q = function(u) qgamma(u, 10, 2), seed = 1,
  estimators = list(
    B_B  = function(b, v) fpa_values(b, method = "bernstein",        bandwidth = "bcv"),
    B_IP = function(b, v) fpa_values(b, method = "indirect_poisson", bandwidth = "bcv"),
    Conv = function(b, v) fpa_values(b, method = "gpv")),
  score = function(fit, b, v) fpa_evaluate_values(fit, pdf, cdf, v))
aggregate(cbind(mise, loglik, ks) ~ estimator, mc, mean)
#>   estimator        mise    loglik         ks
#> 1       B_B 0.007730797 -172.2416 0.09080606
#> 2      B_IP 0.004222934 -192.5422 0.07626960
#> 3      Conv 0.005219661 -179.1400 0.07857900

References

Doosti, H. (2026a). A copula-quantile density approach to affiliated private value auctions. Economics Letters, 261, 112837.

Doosti, H. (2026b). Estimating seller’s expected revenue in first-price sealed-bid auctions via quantile density functions. SSRN preprint. https://doi.org/10.2139/ssrn.7071380

Doosti, H., Dewan, I. and Sbaï, E. (2026). Adaptive wavelet estimation of private-value densities in first-price auctions. SSRN preprint. https://doi.org/10.2139/ssrn.7618223

Doosti, H., Dewan, I. and Talebian, M. (2025). Nonparametric estimation of private value distributions in first-price auctions: Evaluating quantile density function approaches. Economics Letters, 257, 112670.

Guerre, E., Perrigne, I. and Vuong, Q. (2000). Optimal nonparametric estimation of first-price auctions. Econometrica, 68, 525–574.

Hickman, B. R. and Hubbard, T. P. (2015). Replacing sample trimming with boundary correction in nonparametric estimation of first-price auctions. Journal of Applied Econometrics, 30, 739–762.

Ma, J., Marmer, V. and Shneyerov, A. (2019). Inference for first-price auctions with Guerre, Perrigne, and Vuong’s estimator. Journal of Econometrics, 211, 507–538.

Lu, J. and Perrigne, I. (2008). Estimating risk aversion from ascending and sealed-bid auctions: the case of timber auction data. Journal of Applied Econometrics, 23, 871–896.

Zincenko, F. (2024). Estimation and inference of seller’s expected revenue in first-price auctions. Journal of Econometrics, 241, 105734.