En exécutant une fonction (whittST, elle me revoit un résultat "NULL" au lieu d'un nombre calculé préalablement. Ci dessous le code de la fonction, ainsi que ceux de deux autres méthodes appelées. Merci d'avance pour votre aide.
Code : Tout sélectionner
whittST <- function(fBm, Hprel)
{
if (missing(fBm)) fBm <- circFBM(n=500, H=0.6)
if (missing(Hprel))
{
Db4 <- c(0.4829629, -0.8365163, 0.22414386, 0.12940952)
Hprel <- VaPkolST(fBm=fBm, k=2, a=Db4, M=5, 0)$Hols
}
fGn <- c(fBm[1], diff(fBm))
n <- length(fGn)
nstar <- trunc((n - 1)/2)
perio <- (Mod(fft(fGn, inverse = F)))^2/(2 * pi * n)
perio <- perio[1:nstar]
global <- vector("list", 2)
global[[1]] <- perio
global[[2]] <- n
assign ("global", global) ##
## ---------------------
## Criterion to minimize
## ---------------------
Q <- function(Htry)
{
perio <- global[[1]]
n <- global[[2]]
spd <- spdFGN(Htry, n)
result <- (4 * pi)/n * sum(perio/spd)
drop(result)
}
epsilon <- 1.0000000000000001e-05
Hest <- nlminb(objective = Q, start = Hprel, lower = epsilon, upper = 1 - epsilon)$parameters
drop(Hest)
}
whittST(fBm=circFBM(n=1000, H=0.9), Hprel=0.9)
Code : Tout sélectionner
spdFGN <- function(Htry, n)
{
if (missing(Htry)) Htry <- 0.6
if (missing(n)) n <- 500
alpha <- 2 * Htry + 1
nstar <- trunc((n - 1)/2)
clambda <- (sin(pi * Htry) * gamma(alpha))/pi
## -----------------------------------------
## Computation of spd at Fourier frequencies
## -----------------------------------------
j <- 2 * pi * ((-300):300)
spd <- rep(0, nstar)
for (k in (1:nstar))
{
lambda <- (2 * pi * k)/n
stocksum <- sum(abs(j + lambda)^( - alpha))
spd[k] <- clambda * (1 - cos(lambda)) * stocksum
}
## ----------------------
## Renormalization of spd
## ----------------------
theta <- exp(2/n * sum(log(spd)))
spd <- spd/theta
drop(spd)
}
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)
}
}