/
radf_mc.R
166 lines (145 loc) · 4.5 KB
/
radf_mc.R
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
#
#' @importFrom rlang is_scalar_atomic
#' @importFrom doRNG `%dorng%`
#' @importFrom doSNOW registerDoSNOW
radf_mc_ <- function(n, minw, nrep, seed = NULL) {
if(!is.null(dim(n))) {
message_glue("Did you use `data` instead of `n`? Using `NROW(n)` instead.")
n <- NROW(n)
}
assert_n(n)
assert_positive_int(n, greater_than = 5)
assert_positive_int(nrep)
minw <- minw %||% psy_minw(n)
assert_positive_int(minw, greater_than = 2)
show_pb <- getOption("exuber.show_progress")
pb <- set_pb(nrep)
pb_opts <- set_pb_opts(pb)
do_par <- getOption("exuber.parallel")
if (do_par) {
cl <- parallel::makeCluster(getOption("exuber.ncores"), type = "PSOCK")
registerDoSNOW(cl)
on.exit(parallel::stopCluster(cl))
}
`%fun%` <- if (do_par) `%dorng%` else `%do%`#dorng to seed in parallel
set_rng(seed)
results <- foreach(
i = 1:nrep,
.export = c("rls_gsadf", "unroot"),
.combine = "cbind",
.options.snow = pb_opts,
.inorder = FALSE
) %fun% {
if (show_pb && !do_par)
pb$tick()
y <- cumsum(rnorm(n))
yxmat <- unroot(y)
rls_gsadf(yxmat, min_win = minw)
}
n_minw <- n - minw
adf_crit <- results[n_minw + 1, ]
sadf_crit <- results[n_minw + 2, ]
gsadf_crit <- results[n_minw + 3, ]
badf_crit <- results[1:n_minw, ]
bsadf_crit <- results[-c(1:(n_minw + 3)), ]
list(
adf = adf_crit,
sadf = sadf_crit,
gsadf = gsadf_crit,
badf = badf_crit,
bsadf = bsadf_crit) %>%
add_attr(
index = 1:n,
method = "Monte Carlo",
n = n,
minw = minw,
iter = nrep,
seed = get_rng_state(seed),
parallel = do_par
)
}
#' Monte Carlo Critical Values
#'
#' \code{radf_mc_cv} computes Monte Carlo critical values for the recursive unit
#' root tests. \code{radf_mc_distr} computes the distribution.
#'
#' @inheritParams radf
#' @param n A positive integer. The sample size.
#' @param nrep A positive integer. The number of Monte Carlo simulations.
#' @param seed An object specifying if and how the random number generator (rng)
#' should be initialized. Either NULL or an integer will be used in a call to
#' `set.seed` before simulation. If set, the value is saved as "seed" attribute
#' of the returned value. The default, NULL, will not change rng state, and
#' return .Random.seed as the "seed" attribute. Results are different between
#' the parallel and non-parallel option, even if they have the same seed.
#'
#' @return For \code{radf_mc_cv} a list that contains the critical values for ADF,
#' BADF, BSADF and GSADF test statistics. For \code{radf_mc_distr} a list that
#' contains the ADF, SADF and GSADF distributions.
#'
#' @seealso \code{\link{radf_wb_cv}} for wild bootstrap critical values and
#' \code{\link{radf_sb_cv}} for sieve bootstrap critical values
#'
#' @importFrom parallel detectCores makeCluster stopCluster
#' @importFrom foreach foreach %do%
#' @importFrom stats quantile rnorm runif
#' @importFrom lubridate is.Date
#' @importFrom purrr detect_index
#' @export
#'
#' @examples
#' \donttest{
#' # Default minimum window
#' mc <- radf_mc_cv(n = 100)
#'
#' tidy(mc)
#'
#' # Change the minimum window and the number of simulations
#' mc2 <- radf_mc_cv(n = 100, nrep = 600, minw = 20)
#'
#' tidy(mc2)
#'
#' mdist <- radf_mc_distr(n = 100, nrep = 1000)
#'
#' autoplot(mdist)
#' }
radf_mc_cv <- function(n, minw = NULL, nrep = 1000L, seed = NULL) {
pcnt <- c(0.9, 0.95, 0.99)
results <- radf_mc_(n, minw = minw, nrep = nrep, seed = seed)
adf_crit <- quantile(results$adf, probs = pcnt, drop = FALSE)
sadf_crit <- quantile(results$sadf, probs = pcnt, drop = FALSE)
gsadf_crit <- quantile(results$gsadf, probs = pcnt, drop = FALSE)
bsadf_crit <- apply(results$badf, 2, cummax) %>%
apply(1, quantile, probs = pcnt) %>% t()
# values taken from PWY
asy_adf_crit <- rep(
c(-0.44, -0.08, 0.6),
each = nrow(bsadf_crit)
)
badf_crit <- matrix(
asy_adf_crit, ncol = 3,
dimnames = list(NULL, paste0(pcnt*100, "%"))
)
list(
adf_cv = adf_crit,
sadf_cv = sadf_crit,
gsadf_cv = gsadf_crit,
badf_cv = badf_crit,
bsadf_cv = bsadf_crit
) %>%
inherit_attrs(results) %>%
add_class("radf_cv", "mc_cv")
}
#' @rdname radf_mc_cv
#' @inheritParams radf_mc_cv
#' @export
radf_mc_distr <- function(n, minw = NULL, nrep = 1000L, seed = NULL) {
results <- radf_mc_(n, minw = minw, nrep = nrep, seed = seed)
list(
adf_distr = results$adf,
sadf_distr = results$sadf,
gsadf_distr = results$gsadf
) %>%
inherit_attrs(results) %>%
add_class("radf_distr", "mc_distr")
}