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
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
|
#' Wavelet-based Maximum Likelihood Estimation for Seasonal Persistent
#' Processes
#'
#' Parameter estimation for a seasonal persistent (seasonal long-memory)
#' process is performed via maximum likelihood on the wavelet coefficients.
#'
#' The variance-covariance matrix of the original time series is approximated
#' by its wavelet-based equivalent. A Whittle-type likelihood is then
#' constructed where the sums of squared wavelet coefficients are compared to
#' bandpass filtered version of the true spectral density function.
#' Minimization occurs for the fractional difference parameter \eqn{d} and the
#' Gegenbauer frequency \eqn{f_G}, while the innovations variance is
#' subsequently estimated.
#'
#' @usage spp.mle(y, wf, J = log(length(y), 2) - 1, p = 0.01, frac = 1)
#' @usage spp2.mle(y, wf, J = log(length(y), 2) - 1, p = 0.01, dyadic = TRUE, frac = 1)
#' @aliases spp.mle spp2.mle
#' @param y Not necessarily dyadic length time series.
#' @param wf Name of the wavelet filter to use in the decomposition. See
#' \code{\link{wave.filter}} for those wavelet filters available.
#' @param J Depth of the discrete wavelet packet transform.
#' @param p Level of significance for the white noise testing procedure.
#' @param dyadic Logical parameter indicating whether or not the original time
#' series is dyadic in length.
#' @param frac Fraction of the time series that should be used in constructing
#' the likelihood function.
#' @return List containing the maximum likelihood estimates (MLEs) of
#' \eqn{\delta}, \eqn{f_G} and \eqn{\sigma^2}, along with the value of the
#' likelihood for those estimates.
#' @author B. Whitcher
#' @seealso \code{\link{fdp.mle}}
#' @references Whitcher, B. (2004) Wavelet-based estimation for seasonal
#' long-memory processes, \emph{Technometrics}, \bold{46}, No. 2, 225-238.
#' @keywords ts
#' @export spp.mle
spp.mle <- function(y, wf, J=log(length(y),2)-1, p=0.01, frac=1)
{
sppLL <- function(x, y) {
delta <- x[1]
fG <- x[2]
## cat("Parameters are: d =", delta, ", and f =", fG, fill=TRUE)
y.dwpt <- y[[1]]
y.basis <- y[[2]]
n <- y[[3]]
J <- y[[4]]
## Establish the limits of integration for the band-pass variances
a <- unlist(apply(matrix(2^(1:J)-1), 1, seq, from=0, by=1)) /
2^(rep(1:J, 2^(1:J))) / 2
b <- unlist(apply(matrix(2^(1:J)), 1, seq, from=1, by=1)) /
2^(rep(1:J, 2^(1:J))) / 2
## Define some useful parameters for the wavelet packet tree
# n <- length(y)
length.jn <- n / rep(2^(1:J), 2^(1:J))
scale.jn <- rep(2^(1:J+1), 2^(1:J))
## Initialize various parameters for the reduced LL
Basis <- (1:length(y.basis))[y.basis]
bp.var <- numeric(length(Basis))
delta.n <- 100
## Compute the band-pass variances according to \delta and f_G
omega.diag <- NULL
for(i in 1:sum(y.basis)) {
jn <- Basis[i]
bp.var[i] <- bandpass.spp(a[jn], b[jn], delta, fG)
omega.diag <- c(omega.diag,
scale.jn[jn] * rep(bp.var[i], length.jn[jn]))
}
## Compute reduced log-likelihood
rLL <- n * log(1/n * sum(y.dwpt^2 / omega.diag, na.rm=TRUE)) +
sum(length.jn[y.basis] * log(scale.jn[y.basis] * bp.var))
rLL
}
n <- length(y)
x0 <- numeric(2)
## Perform discrete wavelet packet transform (DWPT) on Y
y.dwpt <- dwpt(y, wf, n.levels=J)
n <- length(y)
if(frac < 1) {
for(i in 1:length(y.dwpt)) {
vec <- y.dwpt[[i]]
ni <- length(vec)
j <- rep(1:J, 2^(1:J))[i]
vec[trunc(frac * n/2^j):ni] <- NA
y.dwpt[[i]] <- vec
}
}
y.basis <- as.logical(ortho.basis(portmanteau.test(y.dwpt, p)))
y.dwpt <- as.matrix(unlist(y.dwpt[y.basis]))
## Compute initial estimate of the Gegenbauer frequency
y.per <- per(y - mean(y))
x0[2] <- (0:(n/2)/n)[max(y.per) == y.per]
## Compute initial estimate of the fractional difference parameter
muJ <- (unlist(apply(matrix(2^(1:J)-1), 1, seq, from=0, by=1)) /
2^(rep(1:J, 2^(1:J))) +
unlist(apply(matrix(2^(1:J)), 1, seq, from=1, by=1)) /
2^(rep(1:J, 2^(1:J)))) / 4
y.modwpt <- modwpt(y, wf=wf, n.levels=J)
y.varJ <- rep(2^(1:J), 2^(1:J)) *
unlist(lapply(y.modwpt,
FUN=function(x)sum(x*x,na.rm=TRUE)/length(x[!is.na(x)])))
x0[1] <- min(-0.5 * lsfit(log(abs(muJ[y.basis] - x0[2])),
log(y.varJ[y.basis]))$coef[2], 0.49)
cat(paste("Initial parameters are: delta =", round(x0[1],4),
"freqG =", round(x0[2],4), "\n"))
result <- optim(par=x0, fn=sppLL, method="L-BFGS-B",
lower=c(0.001,0.001), upper=c(0.499,0.499),
control=list(trace=0, fnscale=2),
y=list(y.dwpt, y.basis, n, J))
return(result)
}
spp2.mle <- function(y, wf, J=log(length(y),2)-1, p=0.01,
dyadic=TRUE, frac=1)
{
spp2LL <- function(x, y) {
d1 <- x[1]
f1 <- x[2]
d2 <- x[3]
f2 <- x[4]
## cat("Parameters are: d1 =", round(d1,6), ", and f1 =", round(f1,6),
## ", d2 =", round(d2,6), ", and f2 =", round(f2,6), fill=TRUE)
y.dwpt <- y[[1]]
y.basis <- y[[2]]
n <- y[[3]]
J <- y[[4]]
## Establish the limits of integration for the band-pass variances
a <- unlist(apply(matrix(2^(1:J)-1), 1, seq, from=0, by=1)) /
2^(rep(1:J, 2^(1:J))) / 2
b <- unlist(apply(matrix(2^(1:J)), 1, seq, from=1, by=1)) /
2^(rep(1:J, 2^(1:J))) / 2
## Define some useful parameters for the wavelet packet tree
length.jn <- n / rep(2^(1:J), 2^(1:J))
scale.jn <- rep(2^(1:J+1), 2^(1:J))
## Initialize various parameters for the reduced LL
Basis <- (1:length(y.basis))[y.basis]
bp.var <- numeric(length(Basis))
delta.n <- 100
## Compute the band-pass variances according to \delta and f_G
omega.diag <- NULL
for(i in 1:sum(y.basis)) {
jn <- Basis[i]
bp.var[i] <- bandpass.spp2(a[jn], b[jn], d1, f1, d2, f2)
omega.diag <- c(omega.diag,
scale.jn[jn] * rep(bp.var[i], length.jn[jn]))
}
## Compute reduced log-likelihood
n * log(1/n * sum(y.dwpt^2 / omega.diag, na.rm=TRUE)) +
sum(length.jn[y.basis] * log(scale.jn[y.basis] * bp.var), na.rm=TRUE)
}
n <- length(y)
x0 <- numeric(4)
## Perform discrete wavelet packet transform (DWPT) on Y
y.dwpt <- dwpt(y, wf, n.levels=J)
if(!dyadic) {
for(i in 1:length(y.dwpt)) {
vec <- y.dwpt[[i]]
ni <- length(vec)
j <- rep(1:J, 2^(1:J))[i]
vec[trunc(frac * n/2^j):ni] <- NA
y.dwpt[[i]] <- vec
}
}
y.basis <- as.logical(ortho.basis(portmanteau.test(y.dwpt, p, type="other")))
y.dwpt <- as.vector(unlist(y.dwpt[y.basis]))
## Compute initial estimate of the Gegenbauer frequencies
if(dyadic)
y.per <- per(y - mean(y))
else
y.per <- per(y[1:(frac*n)] - mean(y[1:(frac*n)]))
freq.y <- (0:(frac*n %/% 2))/(frac*n)
x0[2] <- freq.y[max(y.per) == y.per]
x0[4] <- freq.y[max(y.per[freq.y > x0[2] + freq.y[10] |
freq.y < x0[2] - freq.y[10]]) == y.per]
if(x0[2] > x0[4]) {
xx <- x0[2]
x0[2] <- x0[4]
x0[4] <- xx
rm(xx)
}
## Compute initial estimate of the fractional difference parameters
muJ <- (unlist(apply(matrix(2^(1:J)-1), 1, seq, from=0, by=1)) /
2^(rep(1:J, 2^(1:J))) +
unlist(apply(matrix(2^(1:J)), 1, seq, from=1, by=1)) /
2^(rep(1:J, 2^(1:J)))) / 4
y.modwpt <- modwpt(y, wf=wf, n.levels=J)
y.varJ <- rep(2^(1:J), 2^(1:J)) *
unlist(lapply(y.modwpt,
FUN = function(x) sum(x*x,na.rm=TRUE)/length(x[!is.na(x)])))
x0.mid <- (x0[2] + x0[4]) / 2
muJ <- muJ[y.basis]
y.varJ <- y.varJ[y.basis]
x0[1] <- min(-0.5 * lsfit(log(abs(muJ[muJ < x0.mid] - x0[2])),
log(y.varJ[muJ < x0.mid]))$coef[2], 0.49)
x0[3] <- min(-0.5 * lsfit(log(abs(muJ[muJ > x0.mid] - x0[4])),
log(y.varJ[muJ > x0.mid]))$coef[2], 0.49)
cat(paste("Initial parameters: d1 = ", round(x0[1],4),
", f1 = ", round(x0[2],4), ", d2 = ", round(x0[3],4),
", f2 = ", round(x0[4],4), sep=""), fill=TRUE)
result <- optim(par=x0, fn=spp2LL, method="L-BFGS-B",
lower=rep(0.001,4), upper=rep(0.499,4),
control=list(trace=1, fnscale=2),
y=list(y.dwpt, y.basis, n, J))
return(result)
}
|
