sim_sts_quantile_areas.h File Reference

Standardized time series areas of the batched quantile process. More...

#include <algorithm>
#include <cmath>
#include <cstddef>
#include <numeric>
#include <vector>
#include "line/api/sim/sim_types.h"
#include "line/num/number.h"
#include "line/util/error.h"

Include dependency graph for sim_sts_quantile_areas.h:

Go to the source code of this file.

Classes
struct	line::sim::StsQuantileStats< T >
	Batched-quantile statistics of one sample path. More...

Namespaces
namespace	line
namespace	line::sim

Functions
template<class T>
StsQuantileStats< T >	line::sim::sim_sts_quantile_areas (const std::vector< T > &Y, std::size_t b, std::size_t m, double p, double weight=std::sqrt(12.0))
	Standardized time series areas of the batched quantile process.

Detailed Description

Standardized time series areas of the batched quantile process.

Port of matlab/src/api/sim/sim_sts_quantile_areas.m. The b*m observations are split into b nonoverlapping batches of size m; the function returns the signed standardized time series (STS) areas of the quantile-estimation process, the batched quantile estimators, and the three variance-parameter estimators the QUEST procedures build confidence intervals from.

With yhat_p(j,m) the empirical p-quantile of batch j and yhat_p(j,k) that of its first k observations, the STS process of batch j is T_{j,m}(k/m) = (k/sqrt(m)) (yhat_p(j,m) - yhat_p(j,k)), its signed area is A_p(w;j,m) = m^{-1} sum_{k=1}^{m} w(k/m) T_{j,m}(k/m), and the three estimators of sigma_p^2 = lim n Var(ytilde_p(n)) are A_p(w;b,m) = b^{-1} sum_j A_p(w;j,m)^2 (STS area) N_p(b,m) = (b-1)^{-1} m sum_j (yhat_p(j,m)-ytilde_p(n))^2 (NBQ) V_p(w;b,m) = [b A_p(w;b,m) + (b-1) N_p(b,m)] / (2b-1) (combined) where ytilde_p(n) is the full-sample empirical p-quantile over all n = b*m observations. The first two have limiting chi-square laws on b and b-1 degrees of freedom and are asymptotically independent, so the combined estimator carries 2b-1 degrees of freedom and is about sqrt(2) less variable than either component.

WEIGHT FUNCTION. The requirement on w is that int_0^1 w(t)B(t)dt be standard normal for a standard Brownian bridge B; for a constant w = c that variance is c^2/12, so c = sqrt(12) is the normalizing choice and any other constant rescales every area and A_p by (c/sqrt(12))^2. Only constant weights are supported, as in the reference.

COST. The prefix quantiles yhat_p(j,k) are exact order statistics, not a running approximation: a Fenwick tree over the within-batch ranks is advanced one observation at a time and searched by binary lifting, so the whole function is O(b m log m). MATLAB advances the b trees in lockstep to keep its k loop vectorized; here the batches are simply looped, which is the same arithmetic in a different order.

A SINGLE BATCH carries no between-batch degrees of freedom, so N_p and V_p are NaN at b = 1 while areas and bqe stay valid; sim_firquest pools them across replications instead of reading them per replication.

Reference: C. Alexopoulos, D. Goldsman, A. Lolos, K. D. Dingec, J. R. Wilson, "Steady-State Quantile Estimation Using Standardized Time Series", 2020/2023; A. Lolos et al., Proc. Winter Simulation Conference, 2023, theorems 1-3.

Definition in file sim_sts_quantile_areas.h.