|
6 | 6 | #' @param reps numeric; number of replicates to generate. |
7 | 7 | #' @param rho numeric [-1,1] (vectorized); A \code{rho} must be provided, otherwise a blank list will be returned. |
8 | 8 | #' @param type options("spearman", "pearson", "NNScor", "NNSdep"); \code{type = "spearman"}(default) dependence metric desired. |
9 | | -#' @param drift logical; \code{TRUE} default preserves the drift of the original series. |
10 | | -#' @param trim numeric [0,1]; The mean trimming proportion, defaults to \code{trim=0.1}. |
| 9 | +#' @param drift logical; \code{drift = TRUE} (default) preserves the drift of the original series. |
| 10 | +#' @param target_drift numerical; code{NULL} (default) Specifies the desired drift when \code{drift = TRUE}, i.e. a risk-free rate of return. |
| 11 | +#' @param trim numeric [0,1]; The mean trimming proportion, defaults to \code{trim = 0.1}. |
11 | 12 | #' @param xmin numeric; the lower limit for the left tail. |
12 | 13 | #' @param xmax numeric; the upper limit for the right tail. |
13 | 14 | #' @param reachbnd logical; If \code{TRUE} potentially reached bounds (xmin = smallest value - trimmed mean and |
14 | 15 | #' xmax = largest value + trimmed mean) are given when the random draw happens to be equal to 0 and 1, respectively. |
15 | | -#' @param expand.sd logical; If \code{TRUE} the standard deviation in the ensemble is expanded. See \code{expand.sd} in meboot::meboot. |
16 | | -#' @param force.clt logical; If \code{TRUE} the ensemble is forced to satisfy the central limit theorem. See \code{force.clt} in meboot::meboot. |
| 16 | +#' @param expand.sd logical; If \code{TRUE} the standard deviation in the ensemble is expanded. See \code{expand.sd} in \code{meboot::meboot}. |
| 17 | +#' @param force.clt logical; If \code{TRUE} the ensemble is forced to satisfy the central limit theorem. See \code{force.clt} in \code{meboot::meboot}. |
17 | 18 | #' @param scl.adjustment logical; If \code{TRUE} scale adjustment is performed to ensure that the population variance of the transformed series equals the variance of the data. |
18 | 19 | #' @param sym logical; If \code{TRUE} an adjustment is performed to ensure that the ME density is symmetric. |
19 | 20 | #' @param elaps logical; If \code{TRUE} elapsed time during computations is displayed. |
|
59 | 60 | #' boots <- NNS.meboot(AirPassengers, reps=100, rho = 0, xmin = 0) |
60 | 61 | #' |
61 | 62 | #' # Verify correlation of replicates ensemble to original |
62 | | -#' cor(boots["ensemble",], AirPassengers, method = "spearman") |
| 63 | +#' cor(boots["ensemble",]$ensemble, AirPassengers, method = "spearman") |
63 | 64 | #' |
64 | 65 | #' # Plot all replicates |
65 | | -#' matplot(boots["replicates",] , type = 'l') |
| 66 | +#' matplot(boots["replicates",]$replicates , type = 'l') |
66 | 67 | #' |
67 | 68 | #' # Plot ensemble |
68 | | -#' lines(boots["ensemble",], lwd = 3) |
| 69 | +#' lines(boots["ensemble",]$ensemble, lwd = 3) |
69 | 70 | #' } |
70 | 71 | #' @export |
71 | 72 |
|
72 | 73 | NNS.meboot <- function(x, |
73 | | - reps=999, |
74 | | - rho=NULL, |
75 | | - type="spearman", |
76 | | - drift=TRUE, |
77 | | - trim=0.10, |
78 | | - xmin=NULL, |
79 | | - xmax=NULL, |
80 | | - reachbnd=TRUE, |
81 | | - expand.sd=TRUE, force.clt=TRUE, |
82 | | - scl.adjustment = FALSE, sym = FALSE, elaps=FALSE, |
| 74 | + reps = 999, |
| 75 | + rho = NULL, |
| 76 | + type = "spearman", |
| 77 | + drift = TRUE, |
| 78 | + target_drift = NULL, |
| 79 | + trim = 0.10, |
| 80 | + xmin = NULL, |
| 81 | + xmax = NULL, |
| 82 | + reachbnd = TRUE, |
| 83 | + expand.sd = TRUE, force.clt = TRUE, |
| 84 | + scl.adjustment = FALSE, sym = FALSE, elaps = FALSE, |
83 | 85 | digits = 6, |
84 | 86 | colsubj, coldata, coltimes,...) |
85 | 87 | { |
|
206 | 208 | m <- c(matrix2) |
207 | 209 | l <- length(e) |
208 | 210 |
|
209 | | - func <- function(ab, d=drift, ty=type){ |
| 211 | + func <- function(ab, d = drift, ty = type) { |
210 | 212 | a <- ab[1] |
211 | 213 | b <- ab[2] |
212 | | - |
213 | | - if(ty=="spearman" || ty=="pearson"){ |
214 | | - ifelse(d, |
215 | | - (abs(cor((a*m + b*e)/(a + b), e, method = ty) - rho) + |
216 | | - abs(mean((a*m + b*e))/mean(e) - 1) + |
217 | | - abs( cor((a*m + b*e)/(a + b), 1:l) - cor(e, 1:l)) |
218 | | - ), |
219 | | - abs(cor((a*m + b*e)/(a + b), e, method = ty) - rho) + |
220 | | - abs(mean((a*m + b*e))/mean(e) - 1) |
221 | | - ) |
| 214 | + |
| 215 | + # Compute the adjusted ensemble |
| 216 | + combined <- (a * m + b * e) / (a + b) |
| 217 | + |
| 218 | + # Check correlation or dependence structure |
| 219 | + if (ty == "spearman" || ty == "pearson") { |
| 220 | + error <- abs(cor(combined, e, method = ty) - rho) |
| 221 | + } else if (ty == "nnsdep") { |
| 222 | + error <- abs(NNS.dep(combined, e)$Dependence - rho) |
222 | 223 | } else { |
223 | | - if(ty=="nnsdep"){ |
224 | | - ifelse(d, |
225 | | - (abs(NNS.dep((a*m + b*e)/(a + b), e)$Dependence - rho) + |
226 | | - abs(mean((a*m + b*e))/mean(e) - 1) + |
227 | | - abs( NNS.dep((a*m + b*e)/(a + b), 1:l)$Dependence - NNS.dep(e, 1:l)$Dependence) |
228 | | - ), |
229 | | - abs(NNS.dep((a*m + b*e)/(a + b), e)$Dependence - rho) + |
230 | | - abs(mean((a*m + b*e))/mean(e) - 1) |
231 | | - ) |
232 | | - } else { |
233 | | - ifelse(d, |
234 | | - (abs(NNS.dep((a*m + b*e)/(a + b), e)$Correlation - rho) + |
235 | | - abs(mean((a*m + b*e))/mean(e) - 1) + |
236 | | - abs( NNS.dep((a*m + b*e)/(a + b), 1:l)$Correlation - NNS.dep(e, 1:l)$Correlation) |
237 | | - ), |
238 | | - abs(NNS.dep((a*m + b*e)/(a + b), e)$Correlation - rho) + |
239 | | - abs(mean((a*m + b*e))/mean(e) - 1) |
240 | | - ) |
241 | | - } |
| 224 | + error <- abs(NNS.dep(combined, e)$Correlation - rho) |
242 | 225 | } |
243 | 226 |
|
| 227 | + return(error) |
244 | 228 | } |
245 | 229 |
|
246 | | - |
| 230 | + |
247 | 231 | res <- optim(c(.01,.01), func, control=list(abstol = .01)) |
248 | 232 |
|
249 | 233 | ensemble <- (res$par[1]*matrix2 + res$par[2]*ensemble) / (sum(abs(res$par))) |
| 234 | + |
| 235 | + |
| 236 | + # Drift |
| 237 | + orig_coef <- fast_lm(1:n, x)$coef |
| 238 | + orig_intercept <- orig_coef[1] |
| 239 | + orig_drift <- orig_coef[2] |
| 240 | + |
| 241 | + new_coef <- apply(ensemble, 2, function(i) fast_lm(1:n, i)$coef) |
| 242 | + slopes <- new_coef[2,] |
| 243 | + |
| 244 | + if(drift){ |
| 245 | + if(is.null(target_drift)) new_slopes <- (orig_drift - slopes) else new_slopes <- (target_drift - slopes) |
| 246 | + ensemble <- ensemble + t(t(sapply(new_slopes, function(slope) cumsum(rep(slope, n))))) |
| 247 | + |
| 248 | + new_intercepts <- orig_intercept - new_coef[1,] |
| 249 | + ensemble <- sweep(ensemble, 2, new_intercepts, FUN = "+") |
| 250 | + } |
| 251 | + |
| 252 | + |
250 | 253 |
|
251 | 254 | if(identical(ordxx_2, ordxx)){ |
252 | 255 | if(reps>1) ensemble <- t(apply(ensemble, 1, function(x) sample(x, size = reps, replace = TRUE))) |
|
257 | 260 | if(expand.sd) ensemble <- NNS.meboot.expand.sd(x=x, ensemble=ensemble, ...) |
258 | 261 |
|
259 | 262 | if(force.clt && reps > 1) ensemble <- force.clt(x=x, ensemble=ensemble) |
| 263 | + |
| 264 | + |
260 | 265 |
|
261 | 266 | # scale adjustment |
262 | 267 |
|
|
279 | 284 | # Force min / max values |
280 | 285 | if(!is.null(trim[[2]])) ensemble <- apply(ensemble, 2, function(z) pmax(trim[[2]], z)) |
281 | 286 | if(!is.null(trim[[3]])) ensemble <- apply(ensemble, 2, function(z) pmin(trim[[3]], z)) |
| 287 | + |
| 288 | + |
| 289 | + |
282 | 290 |
|
283 | 291 | if(is.ts(x)){ |
284 | 292 | ensemble <- ts(ensemble, frequency=frequency(x), start=start(x)) |
|
287 | 295 | if(reps>1) dimnames(ensemble)[[2]] <- paste("Replicate", 1:reps) |
288 | 296 | } |
289 | 297 |
|
290 | | - |
| 298 | + |
| 299 | + |
291 | 300 | final <- list(x=x, replicates=round(ensemble, digits = digits), ensemble=Rfast::rowmeans(ensemble), xx=xx, z=z, dv=dv, dvtrim=dvtrim, xmin=xmin, |
292 | 301 | xmax=xmax, desintxb=desintxb, ordxx=ordxx, kappa = kappa) |
293 | 302 |
|
|
0 commit comments