CDFs

Cumulative distribution functions (CDFs) represent the probability a variable \(X\) will take a value less than or equal to \(x\). \[F(x) = P(X \leq x)\]

set.seed(123); x = rnorm(100)
ecdf(x)

## Empirical CDF 
## Call: ecdf(x)
##  x[1:100] = -2.3092, -1.9666, -1.6867,  ...,  2.169, 2.1873

P = ecdf(x)
P(0); P(1)

## [1] 0.48

## [1] 0.83

LPM.ratio(degree = 0, target = 0, variable = x); LPM.ratio(degree = 0, target = 1, variable = x)

## [1] 0.48

## [1] 0.83

degree.0.samples = LPM.VaR(percentile = seq(0, 1, length.out = 100), degree = 0, x = x)
degree.0.25.samples = LPM.VaR(percentile = seq(0, 1, length.out = 100), degree = 0.25, x = x)
degree.0.5.samples = LPM.VaR(percentile = seq(0, 1, length.out = 100), degree = 0.5, x = x)
degree.1.samples = LPM.VaR(percentile = seq(0, 1, length.out = 100), degree = 1, x = x)
degree.2.samples = LPM.VaR(percentile = seq(0, 1, length.out = 100), degree = 2, x = x)

head(data.frame(check.names = FALSE, cbind("original x" = sort(x), degree.0.samples, 
                                                          degree.0.25.samples, 
                                                          degree.0.5.samples, 
                                                          degree.1.samples, 
                                                          degree.2.samples)), 10)

     original x degree.0.samples degree.0.25.samples degree.0.5.samples
  1:  -2.309169        -2.309169           -2.309097         -2.3090915
  2:  -1.966617        -1.966617           -1.941190         -1.6935509
  3:  -1.686693        -1.686693           -1.599486         -1.4541494
  4:  -1.548753        -1.548753           -1.382553         -1.2462731
  5:  -1.265396        -1.265396           -1.250823         -1.1453748
  6:  -1.265061        -1.265061           -1.176436         -1.0745440
  7:  -1.220718        -1.220718           -1.119655         -1.0252742
  8:  -1.138137        -1.138137           -1.067793         -0.9868693
  9:  -1.123109        -1.123109           -1.026429         -0.9322105
 10:  -1.071791        -1.071791           -1.014276         -0.8710942
     degree.1.samples degree.2.samples
  1:       -2.3091021       -2.3091170
  2:       -1.4744653       -1.1614908
  3:       -1.2159961       -0.9709972
  4:       -1.0823023       -0.8610192
  5:       -0.9968028       -0.7810300
  6:       -0.9290505       -0.7169770
  7:       -0.8666886       -0.6631888
  8:       -0.8090433       -0.6170691
  9:       -0.7556644       -0.5765608
 10:       -0.7069835       -0.5403318

Bootstrapping (NNS.meboot)

NNS.meboot is based on the maximum entropy bootstrap, available in the R-package meboot. Its rank-based scenario construction is useful for time-series exploration. The default force.clt constrains replicate means using an iid standard-error scale, so preserving rank structure does not by itself provide calibrated sampling uncertainty for serial dependence.

Specified rank-correlation targets explore a chosen range of scenario shapes. They do not enumerate every possible future path. A zero rank-correlation target does not imply independence, and achieved marginal distributions and moments must be checked for the chosen settings.

We will generate 1 replicate of \(X\) for each value of a sequence of \(\rho\) values (the \(ensemble\)), and then plot the results compared to our original \(X\) (black line). NNS.MC is a streamlined wrapper function for this functionality of NNS.meboot.

boots = NNS.MC(x, reps = 1, lower_rho = -1, upper_rho = 1, by = .5)$replicates
reps = do.call(cbind, boots)


matplot(reps, type = "l", col = rainbow(length(boots)))
lines(x, type = "l", lwd = 3, ylim = c(min(reps), max(reps)))

Checking our replicate correlations:

sapply(boots, function(r) cor(r, x, method = "spearman"))

    rho = 1   rho = 0.5     rho = 0  rho = -0.5    rho = -1 
 0.99732373  0.51147915  0.01036904 -0.48720072 -0.98294629 

More replicates and ensembles thereof can be generated for any number of \(\rho\) values.

Note — where the \(\rho\) target lives. The dependence alignment is applied to each replicate individually and in residual space: each replicate’s residuals are mixed between the aligned and anti-aligned orderings of the original residuals so that their dependence matches the requested \(\rho\) (book Chapter 20). For a series without trend, as here, that is also the dependence on the series itself, which the per-replicate correlations above confirm. For a trending series the replicates keep the trend, so their correlation with the original series can stay high even at \(\rho = -1\): on AirPassengers, \(\rho = -1\) replicates have a mean Spearman correlation of about 0.8 with the series. Two things to keep in mind when checking a result: (1) measure it in the same metric that type targets — a type = "NNSdep" target with NNS.dep(...)$Dependence (unsigned), a "spearman"/"pearson" target with rank/linear correlation; a signed correlation taken on an unsigned "NNSdep" target is not comparable to \(\rho\) and can even read negative while the dependence target is met. (2) The ensemble is the per-observation mean of the replicates — a central summary, not a single series carrying the target dependence — so a rank or linear correlation taken directly on the ensemble tends to read higher than the per-replicate \(\rho\) (averaging amplifies the shared order). Calibrate \(\rho\) on the replicates, not the ensemble.

boots = NNS.MC(x, reps = 1, lower_rho = -1, upper_rho = 1, by = .5, target_drift = 0.05)$replicates
reps = do.call(cbind, boots)

plot(x, type = "l", lwd = 3, ylim = c(min(c(x, reps)), max(c(x, reps))))
matplot(reps, type = "l", col = rainbow(length(boots)), add = TRUE)

Simulating a Multivariate Dependence Structure

A marginal rank transformation can transfer the empirical rank arrangement of continuous observations to chosen new marginal distributions. It does not identify a unique copula when marginals have atoms, and a positive-degree partial-moment transform generally changes the original marginal probability law.

Use degree-zero probabilities for this frequency construction. For these continuous samples, subtracting half an observation from the empirical CDF ranks keeps probability levels away from zero and one. Inverting with degree-zero LPM.VaR applies the target empirical quantile function (type7 interpolation).

set.seed(123)
x <- rnorm(1000); y <- rnorm(1000); z <- rnorm(1000)
original.data <- cbind(x, y, z, x_copy = x)
n <- nrow(original.data)
dep.structure <- apply(original.data, 2, function(v) {
  LPM.ratio(degree = 0, target = v, variable = v) - .5 / n
})
new.data <- replicate(ncol(original.data), rnorm(2000, mean = 10, sd = 20))
new.dep.data <- vapply(seq_len(ncol(original.data)), function(j) {
  LPM.VaR(percentile = dep.structure[, j], degree = 0, x = new.data[, j])
}, numeric(n))

cor(original.data, method = "spearman")
cor(new.dep.data, method = "spearman")

The output has one row per original rank tuple, here 1000, even though the new marginal reference samples have 2000 rows. It approximates those reference marginals at the supplied ranks. For tied marginals, specify a tie policy or randomized probability transform within each CDF jump; ordinary inclusive CDF values are not uniform at atoms. Positive-degree transformations remain useful for directional severity modelling, with their changed target stated explicitly.

Alternative Using NNS.meboot

Alternatively, if we wish to keep the simulated values close to the original data, we can apply the NNS.meboot procedure to each of the variables.

We will generate 1 replicate (for brevity) of \(\rho = 0.95\) to our original.data, use their ensemble and note the multivariate dependence among our new.boot.dep.data.

# Apply bootstrap to each variable
new.boot.dep.data = apply(original.data, 2, function(r) NNS.meboot(r, reps = 100, rho = .95))

# Reformat into vectors
boot.ensemble.vectors = lapply(new.boot.dep.data, function(z) unlist(z["ensemble",]))

# Create matrix from vectors
new.boot.dep.matrix = do.call(cbind, boot.ensemble.vectors)

Checking ensemble correlations with original.data:

for(i in 1:4) print(cor(new.boot.dep.matrix[,i], original.data[,i], method = "spearman"))

[1] 0.9885851
[1] 0.9965075
[1] 0.9906302
[1] 0.9901535

Reading these numbers. Here ensemble is the per-observation mean of the replicates, a central summary. Its correlation with the input is not the requested per-replicate target. Alignment is applied to individual residual paths; averaging and trend reconstruction can change the measured correlation. Treat the ensemble as a representative path and inspect achieved correlations for the requested output.

NNS.copula(original.data)
NNS.copula(new.boot.dep.matrix)

[1] 0.4743531
[1] 0.4704039

head(original.data)
head(new.boot.dep.matrix)

               x           y          z           x
[1,] -0.56047565 -0.99579872 -0.5116037 -0.56047565
[2,] -0.23017749 -1.03995504  0.2369379 -0.23017749
[3,]  1.55870831 -0.01798024 -0.5415892  1.55870831
[4,]  0.07050839 -0.13217513  1.2192276  0.07050839
[5,]  0.12928774 -2.54934277  0.1741359  0.12928774
[6,]  1.71506499  1.04057346 -0.6152683  1.71506499
                   x          y          z          x
ensemble1 -0.4268047 -0.7794553 -0.6364458 -0.4642642
ensemble2 -0.2965744 -1.0682197  0.3297265 -0.2531178
ensemble3  1.3302149  0.3054734 -0.4014515  1.4914884
ensemble4  0.2257378  0.3108846  1.0603892  0.1728540
ensemble5  0.4716743 -3.3344967 -0.1917697  0.4309379
ensemble6  1.3984978  1.1881374 -0.5295386  1.5326055