-
Notifications
You must be signed in to change notification settings - Fork 4
/
Copy pathmarglike_NoBreak.R
52 lines (52 loc) · 1.87 KB
/
marglike_NoBreak.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
function (rhos, data, sigma = NULL, k = NULL, type = c("drift",
"both"), norm = TRUE, ...)
{
require(Bolstad)
require(MASS)
data = as.matrix(data)
if (ncol(data) > 1) {
stop("Sorry, 'marglike_NoBreak.R' was written for vector arguments 'data' only not for matrices.\n")
}
like = loglike = lagdiff.star = NULL
y0 = 0
cons = qq = 1
TT = nrow(matrix(data, ncol = 1))
(type = match.arg(arg = type, choices = c("drift", "both"),
several.ok = TRUE))
if (!is.null(k) && k > 1) {
(d1_y = diff(data))
(XX = embed(d1_y, dimension = k))
(xt = embed(cbind(1, 2:TT), dimension = k))
(lagdiff = as.matrix(XX[, -1]))
(lagdiff.star = lagdiff - rep(mean(lagdiff[, 1]), nrow(lagdiff)))
}
else {
xt = cbind(1, 2:TT)
}
(y = matrix(data[-c(1:k)], ncol = 1))
ylag = matrix(data[k:(TT - 1)], ncol = 1)
(xrho = switch(type[1], drift = cbind(xt[, (ncol(xt) - 1)],
lagdiff.star), both = cbind(xt[, (ncol(xt) - 1):ncol(xt)],
lagdiff.star)))
(m.star = t(xrho) %*% xrho)
(inv.m.star = ginv(m.star))
(det.m.star = det(m.star))
if ((-0.01 < det.m.star) && (det.m.star < 0.01)) {
warning("det.m.star==0!! det.m.star<-0.1 has been assigned instead.")
det.m.star = runif(1, 0, 0.1)
}
for (ii in 1:length(rhos)) {
(yrho = y - rhos[ii] * ylag)
(beta.star = inv.m.star %*% (t(xrho) %*% yrho))
(s.star = t(yrho) %*% yrho + (y0^2)/qq - t(beta.star) %*%
m.star %*% beta.star)
(loglikekern = (-(TT - k + 1)/2) * log(s.star) - 0.5 *
log(det.m.star) - 0.5 * log(qq))
loglike = c(loglike, loglikekern)
}
(like = exp(loglike - max(loglike)))
if (norm) {
(cons = sintegral(rhos, like)$value)
}
return(like = like/cons)
}