Pour un algo, j'utilise le package "wavethresh" qui gère notamment les décompositions en ondelettes. j'ai installé le package mais lorsque je fais tourner l'algo, j'ai un message d'erreur qui me dit "could not find function 'wd' ". Cette fonction fait partie de ce package, tout comme "accessD" que j'appelle plus tard. Pouvez-vous m'aider s'il vous plaît ? Ci-joints l'algo appelant les fonctions du package "wavethresh", ainsi que celui d'une fonction appelée pour simuler des trajectoires. Merci d'avance pour votre aide.
Code : Tout sélectionner
waveST <- function(fBm, j1, j2,llplot)
{
n <- length(fBm)
if (missing(fBm)) fBm <- circFBM(n=512, H=0.6)
if ((trunc(log(n)/log(2)) - log(n)/log(2)) != 0)
{
##
## ----------------
## Data duplication
## ----------------
m <- 2
repeat
{
m <- 2 * m
if (m > n) break
}
fBm <- c(fBm, fBm[(n - 1) : (2 * n - m)])
}
if (missing(j1)) j1 <- 2
if (missing(j2)) j2 <- log(length(fBm)) / log(2) - 3
if (missing(llplot)) llplot <- 1
n <- length(fBm)
J <- log(n) / log(2) ##
## ----------------------------------------------
## Wavelet decomposition of fGn and log-scalogram
## ----------------------------------------------
fGnwd <- wd(c(fBm[1], diff(fBm)), 3, family = "DaubExPhase")
scalog <- rep(0, j2 - j1 + 1)
for(j in (j1:j2))
{
## empirical variance of wavelet coefficients
scalog[j - j1 + 1] <- mean(accessD(fGnwd, j)^2)
}
##
## -----------------------------------------------------------
## epsilon_j = log(scalog[j]) / log(2) - (2H + 1)j - log2(K_H)
## estimations of E(epsilon_j) and Var(epsilon_j)
## -----------------------------------------------------------
moy.eps <- c(-0.83527463, -0.39006796, -0.18780535, 0.092044036,
-0.045553664, -0.022659505, -0.011300406, -0.0056428654,
-0.0028195982, -0.0014093405, -0.00070455559, -0.00035224913,
-0.0001761174, -8.8056909e-05, -4.4028006e-05, -2.2013891e-05)
var.eps <- c(3.4237147, 1.3423458, 0.5907403, 0.27710725, 0.13423536,
0.066069658, 0.032776787, 0.016324379, 0.0081462477,
0.0040691462, 0.0020335796, 0.0010165415, 0.00050820866,
0.00025408878, 0.00012704048, 6.3519251e-05)
resolution <- j1:j2
moy.eps <- moy.eps[resolution]
var.eps <- var.eps[resolution] ##
## -----------------------------------------------------------
## Weighted linear regression of log2-scalogram on log2-scales
## -----------------------------------------------------------
Reg <- lsfit(resolution, log(scalog) / log(2) - moy.eps, 1/var.eps)
Hest <- 0.5 * (1 - Reg$coef[2]) ##
## ------------------------------------------
## Bias brought by non-linearity of logarithm
## ------------------------------------------
A <- resolution - mean(resolution)
normA <- as.vector(t(A) %*% A)
bias <- t(A) %*% moy.eps / normA ##
## ------------
## Log-log-plot
## ------------
if(llplot == 1)
{
Hchar <- as.character(round(Hest, 4))
Hleg <- paste(c("Hest", Hchar), collapse = " = ")
par(mfrow = c(1, 1))
plot(resolution, log(scalog)/log(2) - moy.eps, main = "Regression of log2-scalogram on log2-scales")
abline(Reg)
if (Hest < 1/2)
{
yaxis <- max(log(scalog) / log(2) - moy.eps)
}
else yaxis <- (log(scalog[j2 - 2]) / log(2) - moy.eps)
legend(min(resolution), yaxis, c("", Hleg, ""))
}
drop(list(Hest = Hest, bias = bias))
}
waveST(circFBM(n=10000, H=0.62), j1=2, j2=7)
Code : Tout sélectionner
circFBM <- function(n, H, plotfBm)
{
if(missing(n)) n <- 500
if(missing(H)) H <- 0.6
if(missing(plotfBm)) plotfBm <- 1
## -------------------------------------------------------------------
## First line of the circulant matrix, C, built via covariances of fGn
## -------------------------------------------------------------------
lineC <- function(n, H, m)
{
k <- 0:(m - 1)
H2 <- 2 * H
v <- (abs((k - 1)/n)^H2 - 2 * (k/n)^H2 + ((k + 1)/n)^H2)/2
ind <- c(0:(m/2 - 1), m/2, (m/2 - 1):1)
v <- v[ind + 1]
drop(v)
}
## ---------------------
## Next power of two > n
## ---------------------
m <- 2
repeat
{
m <- 2 * m
if (m >= (n - 1))
break
}
stockm <- m ##
## ---------------------------------------------------------------------
## Research of the power of two (<2^18) such that C is definite positive
## ---------------------------------------------------------------------
repeat
{
m <- 2 * m
eigenvalC <- lineC(n, H, m)
eigenvalC <- Re(fft(c(eigenvalC), inverse = F))
if ((all(eigenvalC > 0)) | (m > 2^17))
break
}
if(m > 2^17)
{
cat("----> exact method, impossible!!", fill = T)
cat("----> can't find m such that C is definite positive", fill = T)
break
}
else
{
## ----------------------------------------------------
## Simulation of W=(Q)^t Z, where Z leads N(0, I_m) and
## (Q)_{jk} = m^(-1/2) exp(-2i pi jk/m)
## ----------------------------------------------------
ar <- rnorm(m/2 + 1)
ai <- rnorm(m/2 + 1)
ar[1] <- sqrt(2) * ar[1]
ar[(m/2 + 1)] <- sqrt(2) * ar[(m/2 + 1)]
ai[1] <- 0
ai[(m/2 + 1)] <- 0
ar <- c(ar[c(1:(m/2 + 1))], ar[c((m/2):2)])
aic <- - ai
ai <- c(ai[c(1:(m/2 + 1))], aic[c((m/2):2)])
W <- complex(real = ar, imaginary = ai) ##
## -------------------------
## Reconstruction of the fGn
## -------------------------
W <- (sqrt(eigenvalC)) * W
fGn <- fft(W, inverse = F)
fGn <- (1/(sqrt(2 * m))) * fGn
fGn <- Re(fGn[c(1:n)])
fBm <- cumsum(fGn)
fBm[1] <- 0 ##
## -----------
## Plot of fBm
## -----------
if(plotfBm == 1)
{
par(mfrow = c(1, 1))
time <- (0:(n - 1)) / n
Nchar <- as.character(n)
Nleg <- paste(c("N=", Nchar), collapse = " ")
Hchar <- as.character (round(H, 3))
Hleg <- paste(c(", H=", Hchar), collapse = " ")
NHleg <- paste(c(Nleg, Hleg), collapse = " ")
leg <- paste(c("Path of a fractional Browian motion (Wood and Chan) ---- parameters", NHleg), collapse = " : ")
plot(time, fBm, type = "l", main = leg)
}
drop(fBm)
}
}
circFBM(n=1000, H=0.36)