/
simulate.R
129 lines (116 loc) · 3.87 KB
/
simulate.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
#' Simulate LTM model
#'
#' Simulate LTM model using many
#'
#' @param ns number of times
#' @param nk number of covariates
#' @param ni number of series
#' @param vmu vector mu
#' @param mPhi phi diagonal matrix with the parameters
#' @param mSigs sigma eta vector
#' @param dsig general sigma
#' @param vd threshold parameter
#' @param alpha intercept
#'
#' @return List containing the generated y, x, beta and thresholded beta.
#'
#' @references
#' Nakajima, Jouchi, and Mike West. "Bayesian analysis of latent threshold
#' dynamic models." Journal of Business & Economic
#' Statistics 31.2 (2013): 151-164.
#'
#' @examples
#' # Generates 10 series, each one with 500 observations and 2 regressors.
#'
#' d_sim <- ltm_sim(
#' ns = 500, nk = 2, ni = 10,
#' vmu = matrix(c(.5,.5), nrow = 2),
#' mPhi = diag(2) * c(.99, .99),
#' mSigs = c(.1,.1),
#' dsig = .15,
#' vd = matrix(c(.4,.4), nrow = 2),
#' alpha = 0
#' )
#'
#' str(d_sim)
#'
#' @export
ltm_sim <- function(ns, nk, ni, vmu, mPhi, mSigs, dsig, vd, alpha) {
# variavel preditora
mx <- array(stats::runif(ns*nk*ni)-.5, dim = c(ni, ns, nk))
# betas variando no tempo
mb <- matrix(0, nrow = ns, ncol = nk)
# valor do primeiro beta
mb[1,] <- t(vmu + (chol(solve(diag(nk) - mPhi^2)) * mSigs) %*% stats::rnorm(nk))
# equação do beta
for (i in seq_len(ns-1)) {
mb[i+1,] <- t(vmu + mPhi %*% (t(mb[i,,drop=F]) - vmu) + mSigs * stats::rnorm(nk))
}
# beta zerado
mb_zerado <- mb * t(apply(abs(mb), 1, function(x) x > vd))
# resposta
xbs <- t(apply(mx, 1, function(z) Rfast::rowsums(z * mb_zerado)))
vy <- alpha + xbs + dsig * stats::rnorm(ns*ni)
# resultados
list(vy = vy, mx = mx, mb = mb, mb_zerado = mb_zerado)
}
ltm_pred <- function(post, newdata) {
n_iter <- nrow(post)
nm <- colnames(post)
ni <- sum(grepl("alpha", nm))
nk <- sum(grepl("phi", nm))
nt <- dim(newdata)[2]
soma <- array(0, dim = c(ni, nt))
for(iter in seq_len(n_iter)) {
# recuperando parametros
post_i <- post[iter, ]
mb <- matrix(post_i[grepl("beta", nm)], ncol = nk)
mPhi <- post_i[grepl("phi", nm)] * diag(nk)
mSigs <- post_i[grepl("sig_eta", nm)]
vd <- post_i[grepl("d\\[", nm)]
dsig <- post_i[grepl("sig\\[", nm)]
alpha <- post_i[grepl("alpha\\[", nm)]
vmu <- post_i[grepl("mu\\[", nm)]
# betas variando no tempo
mb_new <- matrix(0, nrow = nt, ncol = nk)
# valor do primeiro beta
mb_new[1,] <- t(vmu + (chol(solve(diag(nk) - mPhi^2)) * mSigs) %*% stats::rnorm(nk))
# equacao do beta
for (i in seq_len(nt-1)) {
mb_new[1,] <- t(vmu + (chol(solve(diag(nk) - mPhi^2)) * mSigs) %*% stats::rnorm(nk))
}
# beta zerado
mb_zerado <- mb_new * t(apply(abs(mb_new), 1, function(x) x > vd))
# resposta
xbs <- t(apply(newdata, 1, function(z) Rfast::rowsums(z * mb_zerado)))
vy <- alpha + xbs + dsig * stats::rnorm(nt*ni)
soma <- soma + vy
}
soma / n_iter
}
ltm_pred_iter <- function(iter, post, nm, ni, nk, newdata) {
# recuperando parametros
post_i <- post[iter, ]
mb <- matrix(post_i[grepl("beta", nm)], ncol = nk)
mPhi <- post_i[grepl("phi", nm)] * diag(nk)
mSigs <- post_i[grepl("sig_eta", nm)]
vd <- post_i[grepl("d\\[", nm)]
dsig <- post_i[grepl("sig\\[", nm)]
alpha <- post_i[grepl("alpha\\[", nm)]
vmu <- post_i[grepl("mu\\[", nm)]
ns_future <- dim(newdata)[2]
# betas variando no tempo
mb_new <- matrix(0, nrow = ns_future, ncol = nk)
# valor do primeiro beta
mb_new[1,] <- t(vmu + (chol(solve(diag(nk) - mPhi^2)) * mSigs) %*% stats::rnorm(nk))
# equação do beta
for (i in seq_len(ns_future-1)) {
mb_new[i+1,] <- t(vmu + mPhi %*% (t(mb_new[i,,drop=F]) - vmu) + mSigs * stats::rnorm(nk))
}
# beta zerado
mb_zerado <- mb_new * t(apply(abs(mb_new), 1, function(x) x > vd))
# resposta
xbs <- t(apply(newdata, 1, function(z) Rfast::rowsums(z * mb_zerado)))
vy <- alpha + xbs + dsig * stats::rnorm(ns_future*ni)
vy
}