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.

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.

rb <- fpa_revenue_boot(timber$n3, r = seq(30, 100, by = 10), B = 30, v_s = 20,
                       method = "bernstein", bandwidth = "wbcv", m_range = c(10, 200))
rb
#> Expected revenue with bootstrap inference (B = 30 , resampling bids )
#>   pointwise CI and uniform band at level 0.95
#>       r estimate      se ci_lower ci_upper band_lower band_upper
#> 30   30  97.5274 17.0632  72.1777 122.8771    72.0780   122.9769
#> 40   40  97.3537 17.2603  71.7524 122.9549    71.6103   123.0971
#> 50   50  96.7570 17.9727  70.1384 123.3756    69.9511   123.5628
#> 60   60  96.6681 18.4772  69.2425 124.0936    69.1097   124.2265
#> 70   70  92.2286 22.8269  59.4457 125.0115    58.1827   126.2745
#> 80   80  89.7598 23.7212  56.0785 123.4411    54.3800   125.1396
#> 90   90  87.0667 24.9620  51.5811 122.5523    49.8365   124.2970
#> 100 100  87.2740 25.0723  52.2334 122.3147    49.8792   124.6689
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.007929796 -171.8784 0.09126195
#> 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. Working paper.

Doosti, H., Dewan, I. and Sbaï, E. (2026). Adaptive wavelet estimation of private-value densities in first-price auctions. Working paper.

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.