Nothing
R/whittle.R
In fArma: Rmetrics - Modelling ARMA Time Series Processes
## # This program is free software; you can redistribute it and/or modify
## # it under the terms of the GNU General Public License as published by
## # the Free Software Foundation; either version 2, or (at your option)
## # any later version.
## #
## # This program is distributed in the hope that it will be useful, but
## # WITHOUT ANY WARRANTY; without even the implied warranty of
## # MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
## # General Public License for more details.
## #
## # A copy of the GNU General Public License is available via WWW at
## # http://www.gnu.org/copyleft/gpl.html. You can also obtain it by
## # writing to the Free Software Foundation, Inc., 59 Temple Place,
## # Suite 330, Boston, MA 02111-1307 USA.
## ################################################################################
## # FUNCTIONS: WHITTLE ESTIMATOR:
## # whittleFit Whittle Estimator
## # .CetaFGN Internal Functions ...
## # .CetaARIMA
## # .Qeta
## # .fspecFGN
## # .ffourier.FGN.est
## # .FGN.spectrum
## # .FGN.B.est.adjust
## # .FGN.B.est
## # .fspecARIMA
## # .per
## # .Qmin2
## # .whittle
## ################################################################################
## ################################################################################
## # DESCRIPTION:
## # The functions are reimplemented from the appendix of J. Beran "Statistics
## # for long-memory processes", Chapman and Hall 1984
## # LICENSE:
## # Permission is hereby given to StatLib to redistribute this software.
## # The software can be freely used for non-commercial purposes, and can
## # be freely distributed for non-commercial purposes only.
## # AUTHORS:
## # Jan Beran <jberan@iris.rz.uni-konstanz.de>
## # Modified: Martin Maechler <maechler@stat.math.ethz.ch>
## # Modified: Diethelm Wuertz <wuertz@itp.phys.ethz.ch> for this R-Port
## whittleFit =
## function(x, order = c(1, 1), subseries = 1, method = c("fgn", "farma"),
## trace = FALSE, spec = FALSE, title = NULL, description = NULL)
## { # A function implemented by Diethelm Wuertz
## # Description:
## # Minimizes an approximate log-likelihood function applied to the
## # spectral density to obtain an estimate of the parameters of a
## # process.
## # Details:
## # Function and programs for the calculation of Whittle's estimator
## # and the goodness of fit statistic as defined in Beran (1992). The
## # models are fractional Gaussian noise or fractional Arima. The data
## # series may be divided into subseries for which the parameters are
## # fitted separately.
## #
## # There are several options for using the Whittle estimator. Some are
## # described below.
## # 1. One can optionally subdivide the series into "subseries".
## # 3. One can output the periodogram by "spec".
## # 4. One can "trace" intermediate minimization results.
## # 5. The "model" can be either farma or fgn.
## # 6. If the model is farma, the "order" has to be specified.
## # 7. The starting value of H for the minimization procedure is "h".
## # 8. "ar" and "ma" are starting values of the time series coefficients.
## # (Length of vectors should be the same as p and q).
## # FUNCTION:
## # Settings:
## data = list(x = x)
## x = as.vector(x)
## # Start Values:
## h = 0.7
## ar = rep(0.5, length = order[1]) / order[1]
## ma = rep(0.5, length = order[2]) / order[2]
## # Estimate:
## if(trace) cat("Iteration Path:\n")
## result = .whittle(xinput = x, nsub = subseries, model = method[1],
## pp = order[1], qq = order[2], h = h, ar = ar, ma = ma, out = trace,
## spec = spec)[[1]]
## result$H = result$par
## result$par = NULL
## # Add:
## if(is.null(title)) title = "Hurst Exponent from Whittle Estimator"
## if(is.null(description)) description = description()
## # Return Value:
## new("fHURST",
## call = match.call(),
## method = paste(method[1], "whittle"),
## hurst = result,
## parameter = list(subseries = subseries, order = order,
## h = h, ar = ar, ma = ma),
## data = data,
## fit = result,
## plot = list(doplot = FALSE),
## title = title,
## description = description
## )
## }
## ################################################################################
## # Functions to make this function independent from Beran's code
## .CetaFGN =
## function(eta)
## { # A function implemented by Diethelm Wuertz
## # Description:
## # Computes covariance matrix of hat{eta} for fGn
## # Author:
## # Jan Beran; modified: Martin Maechler Sep 95, Diethelm Wuertz 2004
## # FUNCTION:
## # Settings:
## M = length(eta)
## # Size of steps in Riemann sum: 2*pi/m
## m = 10000
## mhalfm = trunc((m-1)/2)
## # Size of delta for numerical calculation of derivative
## delta = 1.0e-9
## # Partial derivatives of log f (at each Fourier frequency)
## lf = matrix(1, ncol = M, nrow = mhalfm)
## f0 = .fspecFGN(eta,m)$fspec
## for (j in (1:M)) {
## etaj = eta
## etaj[j] = etaj[j] + delta
## fj = .fspecFGN(etaj, m)$fspec
## lf[,j] = log(fj/f0)/delta }
## # Calculate D:
## Djl = matrix(1,ncol = M, nrow = M)
## for (j in (1:M)) {
## for(l in (1:M)) {
## Djl[j,l] = 2*2*pi/m*sum(lf[,j]*lf[,l])
## }
## }
## ans = drop(matrix(4*pi*solve(Djl), ncol = M, nrow = M, byrow = TRUE))
## # Return Value:
## ans
## }
## # ------------------------------------------------------------------------------
## .CetaARIMA =
## function(eta, p, q)
## { # A function implemented by Diethelm Wuertz
## # Description:
## # Computes ovariance matrix of hat{eta} for fractional ARIMA
## # Author:
## # Jan Beran; modified: Martin Maechler, Diethelm Wuertz 2004
## # FUNCTION:
## # Settings:
## M = length(eta)
## # Size of steps in Riemann sum: 2*pi/m:
## m = 10000
## mhalfm = trunc((m-1)/2)
## # Size of delta for numerical calculation of derivative:
## delta = 1.0e-9
## # partial derivatives of log f (at each Fourier frequency)
## lf = matrix(1, ncol = M, nrow = mhalfm)
## f0 = .fspecARIMA(eta, p, q, m)$fspec
## for (j in (1:M)) {
## etaj = eta
## etaj[j] = etaj[j]+delta
## fj = .fspecARIMA(etaj, p, q, m)$fspec
## lf[,j] = log(fj/f0)/delta
## }
## # Calculate D:
## Djl = matrix(1,ncol = M, nrow = M)
## for (j in (1:M)) {
## for (l in (1:M)) {
## Djl[j,l] = 2*2*pi/m*sum(lf[,j]*lf[,l])
## }
## }
## ans = drop(matrix(4*pi*solve(Djl),ncol = M, nrow = M, byrow = TRUE))
## # Return Value:
## ans
## }
## # ------------------------------------------------------------------------------
## .Qeta =
## function(eta)
## { # A function implemented by Diethelm Wuertz
## # Description:
## # Calculation of A, B and Tn = A/B**2
## # where A = 2pi/n sum 2*[I(lambda = j)/f(lambda = j)],
## # B = 2pi/n sum 2*[I(lambda = j)/f(lambda = j)]**2 and
## # the sum is taken over all Fourier frequencies
## # lambda = j = 2pi*j/n (j=1,...,(n-1)/2.
## # f is the spectral density of fractional Gaussian
## # noise or fractional ARIMA(p,d,q) with self-similarity parameter H = h.
## # cov(X(t),X(t+k))=integral(exp(iuk)f(u)du)
## # Arguments:
## # h
## # (n, nhalfm = trunc[(n-1)/2] and the
## # nhalfm-dimensional GLOBAL vector `yper' must be defined.)
## # Value:
## # list(n=n,h=h,A=A,B=B,Tn=Tn,z=z,pval=pval, theta1=theta1,fspec=fspec)
## # Tn is the goodness of fit test statistic
## # Tn=A/B**2 defined in Beran (1992),
## # z is the standardized test statistic,
## # pval the corresponding p-value P(w>z).
## # theta1 is the scale parameter such that
## # f=theta1*fspec and integral(log[fspec]) = 0.
## # Note:
## # yper[1] must be the periodogram I(lambda = 1) at
## # the frequency 2pi/n (i.e. not the frequency zero !).
## # Author:
## # Jan Beran; modified: Martin Maechler Sep. 95, Diethelm Wuertz 2004
## # FUNCTION:
## # To Suppress No Visible Bindings Warning/Error:
## if(FALSE) { imodel = n = p = yper = NA }
## # Settings:
## h = eta[1]
## if(imodel == 1) {
## fspec = .fspecFGN(eta, n)
## theta1 = fspec$theta1
## fspec = fspec$fspec
## } else {
## fspec = .fspecARIMA(eta, p, q, n)
## theta1 = fspec$theta1
## fspec = fspec$fspec
## }
## yf = yper/fspec
## yfyf = yf**2
## A = 2*(2*pi/n)*sum(yfyf)
## B = 2*(2*pi/n)*sum(yf)
## Tn = A/(B**2)
## z = sqrt(n)*(pi*Tn-1)/sqrt(2)
## pval = 1-pnorm(z)
## theta1 = B/(2*pi)
## fspec = fspec
## Qresult = list(n = n, h = h, eta = eta, A = A,B = B, Tn = Tn,
## z = z, pval = pval, theta1 = theta1, fspec = fspec)
## ans = drop(Qresult)
## # Return value:
## ans
## }
## # ------------------------------------------------------------------------------
## .fspecFGN =
## function(eta, m)
## { # A function implemented by Diethelm Wuertz
## # Description:
## # Calculation of the spectral density f of normalized fractional
## # Gaussian noise with self-similarity parameter H=h at the
## # Fourier frequencies 2*pi*j/m (j=1,...,(m-1)).
## # Arguments:
## # m = sample size
## # h = self-similarity parameter
## # Value:
## # list(fspec = fspec, theta1 = theta1)
## # Note:
## # 1. cov(X(t),X(t+k)) = integral[exp(iuk)f(u)du]
## # 2. f = theta1*fspec and integral[log(fspec)] = 0.
## # Author:
## # Jan Beran; modified: Martin Maechler, Diethelm Wuertz 2004
## # FUNCTION:
## # Settings:
## # Taqqu: This Implementation is more efficient than Beran's:
## fspec = .ffourier.FGN.est(eta, m)
## logfspec = log(fspec)
## fint = 2/(m)*sum(logfspec)
## theta1 = exp(fint)
## fspec = fspec/theta1
## ans = drop(list(fspec = fspec, theta1 = theta1))
## # Return Value:
## ans
## }
## # ------------------------------------------------------------------------------
## .ffourier.FGN.est =
## function(H, n)
## { # A function implemented by Diethelm Wuertz
## # Description:
## # Internal Function
## # Author:
## # Jan Beran; modified: Martin Maechler, Diethelm Wuertz 2004
## # FUNCTION:
## # Spectrum:
## ans = .FGN.spectrum((2 * pi * (1:((n - 1)/2)))/n, H)/pi/2
## # Return Value:
## ans
## }
## # ------------------------------------------------------------------------------
## .FGN.spectrum =
## function(lambda, H)
## { # A function implemented by Diethelm Wuertz
## # Description:
## # Internal Function
## # Author:
## # Jan Beran; modified: Martin Maechler, Diethelm Wuertz 2004
## # FUNCTION:
## # Settings:
## ans = 2*sin(pi*H)*gamma(2*H+1)*(1-cos(lambda))*(lambda^(-2*H-1) +
## .FGN.B.est.adjust(lambda, H))
## ans
## }
## # ------------------------------------------------------------------------------
## .FGN.B.est.adjust =
## function(lambda, H)
## { # A function implemented by Diethelm Wuertz
## # Description:
## # Internal Function
## # Author:
## # Jan Beran; modified: Martin Maechler, Diethelm Wuertz 2004
## # FUNCTION:
## # Settings:
## B = .FGN.B.est(lambda, H)
## ans = (1.0002-0.000134*lambda) * (B-2^(-7.65*H-7.4))
## ans
## }
## # ------------------------------------------------------------------------------
## .FGN.B.est =
## function(lambda, H)
## { # A function implemented by Diethelm Wuertz
## # Description:
## # Internal Function
## # Author:
## # Jan Beran; modified: Martin Maechler, Diethelm Wuertz 2004
## # FUNCTION:
## # Settings:
## d = -(2*H+1)
## dprime = -2*H
## a = function(lambda, k) { 2 * k * pi+lambda }
## b = function(lambda, k) { 2 * k * pi-lambda }
## a1 = a(lambda, 1); b1 = b(lambda, 1)
## a2 = a(lambda, 2); b2 = b(lambda, 2)
## a3 = a(lambda, 3)
## b3 = b(lambda, 3)
## a4 = a(lambda, 4)
## b4 = b(lambda, 4)
## ans = a1^d+b1^d+a2^d+b2^d+a3^d+b3^d+
## (a3^dprime+b3^dprime+a4^dprime+b4^dprime)/(8*pi*H)
## ans
## }
## # ------------------------------------------------------------------------------
## .fspecARIMA =
## function(eta, p, q, m)
## { # A function implemented by Diethelm Wuertz
## # Description:
## # Internal Function
## # Author:
## # Jan Beran; modified: Martin Maechler, Diethelm Wuertz 2004
## # FUNCTION:
## # Settings:
## h = eta[1]
## phi = c()
## psi = c()
## mhalfm = trunc((m-1)/2)
## x = 2*pi/m*(1:mhalfm)
## # Calculation of f at Fourier frequencies
## far = (1:mhalfm)/(1:mhalfm)
## fma = (1:mhalfm)/(1:mhalfm)
## if(p > 0) {
## phi = cbind(eta[2:(p+1)])
## cosar = cos(cbind(x) %*% rbind(1:p))
## sinar = sin(cbind(x) %*% rbind(1:p))
## Rar = cosar %*% phi
## Iar = sinar %*% phi
## far = (1-Rar)**2 + Iar**2 }
## if(q > 0) {
## psi = cbind(eta[(p+2):(p+q+1)])
## cosar = cos(cbind(x) %*% rbind(1:q))
## sinar = sin(cbind(x) %*% rbind(1:q))
## Rar = cosar %*% psi
## Iar = sinar %*% psi
## fma = (1+Rar)**2 + Iar**2 }
## fspec = fma/far*sqrt((1-cos(x))**2 + sin(x)**2)**(1-2*h)
## theta1 = 1/(2*pi)
## ans = list(fspec = fspec, theta1 = theta1)
## ans
## }
## # ------------------------------------------------------------------------------
## .per =
## function(z)
## { # A function implemented by Diethelm Wuertz
## # Description:
## # Internal Function
## # Author:
## # Jan Beran; modified: Martin Maechler, Diethelm Wuertz 2004
## # FUNCTION:
## # Settings:
## n = length(z)
## ans = (Mod(fft(z))**2/(2*pi*n))[1:(n %/% 2 + 1)]
## # Return Value:
## ans
## }
## # ------------------------------------------------------------------------------
## .Qmin2 =
## function(etatry)
## { # A function implemented by Diethelm Wuertz
## # Description:
## # Internal Function
## # FUNCTION:
## # Compute:
## ans = .Qeta(etatry)$B
## assign("bBb", ans, pos = 1)
## # Return Value:
## ans
## }
## # ------------------------------------------------------------------------------
## # Internal Function: WHITTLE ESTIMATOR:
## .whittle =
## function(xinput, nsub = 1, model = c("farma", "fgn"), pp = 1,
## qq = 1, h = 0.5, ar = c(0.5), ma = c(0.5), out = TRUE, spec = FALSE)
## { # A function implemented by Diethelm Wuertz
## # Description:
## # Internal Function
## # Author:
## # Jan Beran; modified: Martin Maechler, Diethelm Wuertz 2004
## # FUNCTION:
## # To suppress No Visible Bindings Warning/Error:
## if(FALSE) { n = p = q = imodel = out = yper = bBb = NA }
## # Settings:
## model = model[1]
## assign("out", out, pos = 1)
## nmax = length(xinput)
## startend = c(1, nmax)
## istart = startend[1]
## iend = startend[2]
## nloop = nsub
## assign("n", trunc((iend - istart+1)/nloop), pos = 1)
## nhalfm = trunc((n - 1)/2)
## if(model == "farma")
## assign("imodel", 2, pos = 1)
## if(model == "fgn")
## assign("imodel", 1, pos = 1)
## assign("p", 0, pos = 1)
## assign("q", 0, pos = 1)
## p <- q <- 0
## if(imodel == 2) {
## assign("p", pp, pos = 1); assign("q", qq, pos = 1)
## } else {
## assign("p", 0, pos=1); assign("q", 0, pos = 1)
## }
## eta = c(h)
## if(p > 0) eta[2:(p+1)] = ar
## if(q > 0) eta[(p+2):(p+q+1)] = ma
## M = length(eta) #loop
## thetavector = c()
## i0 = istart
## flax = vector("list", nloop)
## # Necessary to make nsub/nloop work. VT.
## for (iloop in (1:nloop)) {
## h = max(0.2, min(h, 0.9))
## eta[1] = h
## i1 = i0+n - 1
## y = xinput[i0:i1]
## # Standardize Data:
## vary = var(y)
## y = (y - mean(y))/sqrt(var(y))
## # Periodogram of the Data:
## if(spec) {
## assign("yper", .per(y)[2:(nhalfm+1)], pos = 1)
## } else {
## assign("yper", .per(y)[2:(nhalfm+1)], pos = 1)
## }
## s = 2*(1-h)
## etatry = eta
## # Modified to make optim not give incorrect result. VT
## if(imodel == 1) {
## result = optim(par = etatry, fn = .Qmin2,
## method = "L-BFGS-B", lower = 0, upper = 0.999)
## } else {
## result = optim(par = etatry, fn = .Qmin2)
## }
## eta = result$par
## sturno = result$message
## theta1 = .Qeta(eta)$theta1
## theta = c(theta1, eta)
## thetavector = c(thetavector, theta)
## # Calculate goodness of fit statistic
## Qresult = .Qeta(eta) #output
## M = length(eta)
## if(imodel == 1) {
## SD = .CetaFGN(eta)
## SD = matrix(SD, ncol = M, nrow = M, byrow = TRUE) / n
## } else {
## # Changed to eliminate crashing in solve.qr in CetaARIMA. VT
## cat("M =", M, "\n")
## if(M > 2) {
## for (i in 3:M) {
## for (j in 2:(i-1)) {
## temp = eta[i]+eta[j]
## if(abs(temp) < 0.0001) {
## cat("Problem with estimating confidence intervals,",
## "parameter ", i, "and parameter ", j,
## "are the same, eliminating.\n")
## eta = eta[ - i]
## eta = eta[ - j]
## M = M - 2
## p = p - 1
## q = q - 1
## }
## }
## }
## }
## SD = .CetaARIMA(eta, p, q)
## SD = matrix(SD, ncol = M, nrow = M, byrow = TRUE)/n
## }
## Hlow = eta[1] - 1.96 * sqrt(SD[1, 1])
## Hup = eta[1]+1.96 * sqrt(SD[1, 1])
## if(out) {
## cat("theta =", theta, fill = TRUE)
## cat("H =", eta[1], fill = TRUE)
## cat("95%-CI for H: [", Hlow, ",", Hup, "]", fill = TRUE)
## }
## # Changing of the signs of the moving average parameters
## # in order to respect the sign of the Splus convention
## if(q > 0) eta[(p+2):(p+q+1)] = -eta[(p+2):(p+q+1)]
## etalow = c()
## etaup = c()
## for (i in (1:length(eta))) {
## etalow = c(etalow, eta[i] - 1.96 * sqrt(SD[i, i]))
## etaup = c(etaup, eta[i]+1.96 * sqrt(SD[i, i]))
## }
## if(out) {
## cat("95%-CI:", fill = TRUE)
## print(cbind(etalow, etaup), fill = TRUE)
## }
## if(spec) {
## cat("Periodogram is in yper", fill = TRUE)
## assign("fest", Qresult$theta1 * Qresult$fspec, pos = 1)
## cat("Spectral density is in fest", fill = TRUE)
## }
## flax[[iloop]] = list()
## flax[[iloop]]$par = eta
## flax[[iloop]]$sigma2 = bBb * var(xinput)
## flax[[iloop]]$conv.type = sturno
## remove("bBb", pos = 1)
## # Next subseries:
## i0 = i0+n
## # Changing of the signs of the moving average parameters:
## if(q > 0) eta[(p+2):(p+q+1)] = -eta[(p+2):(p+q+1)]
## } # end of nloop
## # Return:
## return(flax)
## # Clean up:
## remove("n", pos = 1)
## remove("p", pos = 1)
## remove("q", pos = 1)
## remove("imodel", pos = 1)
## remove("out", pos = 1)
## if(spec == FALSE) remove("yper", pos = 1)
## }
## ################################################################################
Try the fArma package in your browser
Any scripts or data that you put into this service are public.
fArma documentation built on Sept. 9, 2022, 3:02 p.m.
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
## # This program is free software; you can redistribute it and/or modify
## # it under the terms of the GNU General Public License as published by
## # the Free Software Foundation; either version 2, or (at your option)
## # any later version.
## #
## # This program is distributed in the hope that it will be useful, but
## # WITHOUT ANY WARRANTY; without even the implied warranty of
## # MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
## # General Public License for more details.
## #
## # A copy of the GNU General Public License is available via WWW at
## # http://www.gnu.org/copyleft/gpl.html. You can also obtain it by
## # writing to the Free Software Foundation, Inc., 59 Temple Place,
## # Suite 330, Boston, MA 02111-1307 USA.
## ################################################################################
## # FUNCTIONS: WHITTLE ESTIMATOR:
## # whittleFit Whittle Estimator
## # .CetaFGN Internal Functions ...
## # .CetaARIMA
## # .Qeta
## # .fspecFGN
## # .ffourier.FGN.est
## # .FGN.spectrum
## # .FGN.B.est.adjust
## # .FGN.B.est
## # .fspecARIMA
## # .per
## # .Qmin2
## # .whittle
## ################################################################################
## ################################################################################
## # DESCRIPTION:
## # The functions are reimplemented from the appendix of J. Beran "Statistics
## # for long-memory processes", Chapman and Hall 1984
## # LICENSE:
## # Permission is hereby given to StatLib to redistribute this software.
## # The software can be freely used for non-commercial purposes, and can
## # be freely distributed for non-commercial purposes only.
## # AUTHORS:
## # Jan Beran <jberan@iris.rz.uni-konstanz.de>
## # Modified: Martin Maechler <maechler@stat.math.ethz.ch>
## # Modified: Diethelm Wuertz <wuertz@itp.phys.ethz.ch> for this R-Port
## whittleFit =
## function(x, order = c(1, 1), subseries = 1, method = c("fgn", "farma"),
## trace = FALSE, spec = FALSE, title = NULL, description = NULL)
## { # A function implemented by Diethelm Wuertz
## # Description:
## # Minimizes an approximate log-likelihood function applied to the
## # spectral density to obtain an estimate of the parameters of a
## # process.
## # Details:
## # Function and programs for the calculation of Whittle's estimator
## # and the goodness of fit statistic as defined in Beran (1992). The
## # models are fractional Gaussian noise or fractional Arima. The data
## # series may be divided into subseries for which the parameters are
## # fitted separately.
## #
## # There are several options for using the Whittle estimator. Some are
## # described below.
## # 1. One can optionally subdivide the series into "subseries".
## # 3. One can output the periodogram by "spec".
## # 4. One can "trace" intermediate minimization results.
## # 5. The "model" can be either farma or fgn.
## # 6. If the model is farma, the "order" has to be specified.
## # 7. The starting value of H for the minimization procedure is "h".
## # 8. "ar" and "ma" are starting values of the time series coefficients.
## # (Length of vectors should be the same as p and q).
## # FUNCTION:
## # Settings:
## data = list(x = x)
## x = as.vector(x)
## # Start Values:
## h = 0.7
## ar = rep(0.5, length = order[1]) / order[1]
## ma = rep(0.5, length = order[2]) / order[2]
## # Estimate:
## if(trace) cat("Iteration Path:\n")
## result = .whittle(xinput = x, nsub = subseries, model = method[1],
## pp = order[1], qq = order[2], h = h, ar = ar, ma = ma, out = trace,
## spec = spec)[[1]]
## result$H = result$par
## result$par = NULL
## # Add:
## if(is.null(title)) title = "Hurst Exponent from Whittle Estimator"
## if(is.null(description)) description = description()
## # Return Value:
## new("fHURST",
## call = match.call(),
## method = paste(method[1], "whittle"),
## hurst = result,
## parameter = list(subseries = subseries, order = order,
## h = h, ar = ar, ma = ma),
## data = data,
## fit = result,
## plot = list(doplot = FALSE),
## title = title,
## description = description
## )
## }
## ################################################################################
## # Functions to make this function independent from Beran's code
## .CetaFGN =
## function(eta)
## { # A function implemented by Diethelm Wuertz
## # Description:
## # Computes covariance matrix of hat{eta} for fGn
## # Author:
## # Jan Beran; modified: Martin Maechler Sep 95, Diethelm Wuertz 2004
## # FUNCTION:
## # Settings:
## M = length(eta)
## # Size of steps in Riemann sum: 2*pi/m
## m = 10000
## mhalfm = trunc((m-1)/2)
## # Size of delta for numerical calculation of derivative
## delta = 1.0e-9
## # Partial derivatives of log f (at each Fourier frequency)
## lf = matrix(1, ncol = M, nrow = mhalfm)
## f0 = .fspecFGN(eta,m)$fspec
## for (j in (1:M)) {
## etaj = eta
## etaj[j] = etaj[j] + delta
## fj = .fspecFGN(etaj, m)$fspec
## lf[,j] = log(fj/f0)/delta }
## # Calculate D:
## Djl = matrix(1,ncol = M, nrow = M)
## for (j in (1:M)) {
## for(l in (1:M)) {
## Djl[j,l] = 2*2*pi/m*sum(lf[,j]*lf[,l])
## }
## }
## ans = drop(matrix(4*pi*solve(Djl), ncol = M, nrow = M, byrow = TRUE))
## # Return Value:
## ans
## }
## # ------------------------------------------------------------------------------
## .CetaARIMA =
## function(eta, p, q)
## { # A function implemented by Diethelm Wuertz
## # Description:
## # Computes ovariance matrix of hat{eta} for fractional ARIMA
## # Author:
## # Jan Beran; modified: Martin Maechler, Diethelm Wuertz 2004
## # FUNCTION:
## # Settings:
## M = length(eta)
## # Size of steps in Riemann sum: 2*pi/m:
## m = 10000
## mhalfm = trunc((m-1)/2)
## # Size of delta for numerical calculation of derivative:
## delta = 1.0e-9
## # partial derivatives of log f (at each Fourier frequency)
## lf = matrix(1, ncol = M, nrow = mhalfm)
## f0 = .fspecARIMA(eta, p, q, m)$fspec
## for (j in (1:M)) {
## etaj = eta
## etaj[j] = etaj[j]+delta
## fj = .fspecARIMA(etaj, p, q, m)$fspec
## lf[,j] = log(fj/f0)/delta
## }
## # Calculate D:
## Djl = matrix(1,ncol = M, nrow = M)
## for (j in (1:M)) {
## for (l in (1:M)) {
## Djl[j,l] = 2*2*pi/m*sum(lf[,j]*lf[,l])
## }
## }
## ans = drop(matrix(4*pi*solve(Djl),ncol = M, nrow = M, byrow = TRUE))
## # Return Value:
## ans
## }
## # ------------------------------------------------------------------------------
## .Qeta =
## function(eta)
## { # A function implemented by Diethelm Wuertz
## # Description:
## # Calculation of A, B and Tn = A/B**2
## # where A = 2pi/n sum 2*[I(lambda = j)/f(lambda = j)],
## # B = 2pi/n sum 2*[I(lambda = j)/f(lambda = j)]**2 and
## # the sum is taken over all Fourier frequencies
## # lambda = j = 2pi*j/n (j=1,...,(n-1)/2.
## # f is the spectral density of fractional Gaussian
## # noise or fractional ARIMA(p,d,q) with self-similarity parameter H = h.
## # cov(X(t),X(t+k))=integral(exp(iuk)f(u)du)
## # Arguments:
## # h
## # (n, nhalfm = trunc[(n-1)/2] and the
## # nhalfm-dimensional GLOBAL vector `yper' must be defined.)
## # Value:
## # list(n=n,h=h,A=A,B=B,Tn=Tn,z=z,pval=pval, theta1=theta1,fspec=fspec)
## # Tn is the goodness of fit test statistic
## # Tn=A/B**2 defined in Beran (1992),
## # z is the standardized test statistic,
## # pval the corresponding p-value P(w>z).
## # theta1 is the scale parameter such that
## # f=theta1*fspec and integral(log[fspec]) = 0.
## # Note:
## # yper[1] must be the periodogram I(lambda = 1) at
## # the frequency 2pi/n (i.e. not the frequency zero !).
## # Author:
## # Jan Beran; modified: Martin Maechler Sep. 95, Diethelm Wuertz 2004
## # FUNCTION:
## # To Suppress No Visible Bindings Warning/Error:
## if(FALSE) { imodel = n = p = yper = NA }
## # Settings:
## h = eta[1]
## if(imodel == 1) {
## fspec = .fspecFGN(eta, n)
## theta1 = fspec$theta1
## fspec = fspec$fspec
## } else {
## fspec = .fspecARIMA(eta, p, q, n)
## theta1 = fspec$theta1
## fspec = fspec$fspec
## }
## yf = yper/fspec
## yfyf = yf**2
## A = 2*(2*pi/n)*sum(yfyf)
## B = 2*(2*pi/n)*sum(yf)
## Tn = A/(B**2)
## z = sqrt(n)*(pi*Tn-1)/sqrt(2)
## pval = 1-pnorm(z)
## theta1 = B/(2*pi)
## fspec = fspec
## Qresult = list(n = n, h = h, eta = eta, A = A,B = B, Tn = Tn,
## z = z, pval = pval, theta1 = theta1, fspec = fspec)
## ans = drop(Qresult)
## # Return value:
## ans
## }
## # ------------------------------------------------------------------------------
## .fspecFGN =
## function(eta, m)
## { # A function implemented by Diethelm Wuertz
## # Description:
## # Calculation of the spectral density f of normalized fractional
## # Gaussian noise with self-similarity parameter H=h at the
## # Fourier frequencies 2*pi*j/m (j=1,...,(m-1)).
## # Arguments:
## # m = sample size
## # h = self-similarity parameter
## # Value:
## # list(fspec = fspec, theta1 = theta1)
## # Note:
## # 1. cov(X(t),X(t+k)) = integral[exp(iuk)f(u)du]
## # 2. f = theta1*fspec and integral[log(fspec)] = 0.
## # Author:
## # Jan Beran; modified: Martin Maechler, Diethelm Wuertz 2004
## # FUNCTION:
## # Settings:
## # Taqqu: This Implementation is more efficient than Beran's:
## fspec = .ffourier.FGN.est(eta, m)
## logfspec = log(fspec)
## fint = 2/(m)*sum(logfspec)
## theta1 = exp(fint)
## fspec = fspec/theta1
## ans = drop(list(fspec = fspec, theta1 = theta1))
## # Return Value:
## ans
## }
## # ------------------------------------------------------------------------------
## .ffourier.FGN.est =
## function(H, n)
## { # A function implemented by Diethelm Wuertz
## # Description:
## # Internal Function
## # Author:
## # Jan Beran; modified: Martin Maechler, Diethelm Wuertz 2004
## # FUNCTION:
## # Spectrum:
## ans = .FGN.spectrum((2 * pi * (1:((n - 1)/2)))/n, H)/pi/2
## # Return Value:
## ans
## }
## # ------------------------------------------------------------------------------
## .FGN.spectrum =
## function(lambda, H)
## { # A function implemented by Diethelm Wuertz
## # Description:
## # Internal Function
## # Author:
## # Jan Beran; modified: Martin Maechler, Diethelm Wuertz 2004
## # FUNCTION:
## # Settings:
## ans = 2*sin(pi*H)*gamma(2*H+1)*(1-cos(lambda))*(lambda^(-2*H-1) +
## .FGN.B.est.adjust(lambda, H))
## ans
## }
## # ------------------------------------------------------------------------------
## .FGN.B.est.adjust =
## function(lambda, H)
## { # A function implemented by Diethelm Wuertz
## # Description:
## # Internal Function
## # Author:
## # Jan Beran; modified: Martin Maechler, Diethelm Wuertz 2004
## # FUNCTION:
## # Settings:
## B = .FGN.B.est(lambda, H)
## ans = (1.0002-0.000134*lambda) * (B-2^(-7.65*H-7.4))
## ans
## }
## # ------------------------------------------------------------------------------
## .FGN.B.est =
## function(lambda, H)
## { # A function implemented by Diethelm Wuertz
## # Description:
## # Internal Function
## # Author:
## # Jan Beran; modified: Martin Maechler, Diethelm Wuertz 2004
## # FUNCTION:
## # Settings:
## d = -(2*H+1)
## dprime = -2*H
## a = function(lambda, k) { 2 * k * pi+lambda }
## b = function(lambda, k) { 2 * k * pi-lambda }
## a1 = a(lambda, 1); b1 = b(lambda, 1)
## a2 = a(lambda, 2); b2 = b(lambda, 2)
## a3 = a(lambda, 3)
## b3 = b(lambda, 3)
## a4 = a(lambda, 4)
## b4 = b(lambda, 4)
## ans = a1^d+b1^d+a2^d+b2^d+a3^d+b3^d+
## (a3^dprime+b3^dprime+a4^dprime+b4^dprime)/(8*pi*H)
## ans
## }
## # ------------------------------------------------------------------------------
## .fspecARIMA =
## function(eta, p, q, m)
## { # A function implemented by Diethelm Wuertz
## # Description:
## # Internal Function
## # Author:
## # Jan Beran; modified: Martin Maechler, Diethelm Wuertz 2004
## # FUNCTION:
## # Settings:
## h = eta[1]
## phi = c()
## psi = c()
## mhalfm = trunc((m-1)/2)
## x = 2*pi/m*(1:mhalfm)
## # Calculation of f at Fourier frequencies
## far = (1:mhalfm)/(1:mhalfm)
## fma = (1:mhalfm)/(1:mhalfm)
## if(p > 0) {
## phi = cbind(eta[2:(p+1)])
## cosar = cos(cbind(x) %*% rbind(1:p))
## sinar = sin(cbind(x) %*% rbind(1:p))
## Rar = cosar %*% phi
## Iar = sinar %*% phi
## far = (1-Rar)**2 + Iar**2 }
## if(q > 0) {
## psi = cbind(eta[(p+2):(p+q+1)])
## cosar = cos(cbind(x) %*% rbind(1:q))
## sinar = sin(cbind(x) %*% rbind(1:q))
## Rar = cosar %*% psi
## Iar = sinar %*% psi
## fma = (1+Rar)**2 + Iar**2 }
## fspec = fma/far*sqrt((1-cos(x))**2 + sin(x)**2)**(1-2*h)
## theta1 = 1/(2*pi)
## ans = list(fspec = fspec, theta1 = theta1)
## ans
## }
## # ------------------------------------------------------------------------------
## .per =
## function(z)
## { # A function implemented by Diethelm Wuertz
## # Description:
## # Internal Function
## # Author:
## # Jan Beran; modified: Martin Maechler, Diethelm Wuertz 2004
## # FUNCTION:
## # Settings:
## n = length(z)
## ans = (Mod(fft(z))**2/(2*pi*n))[1:(n %/% 2 + 1)]
## # Return Value:
## ans
## }
## # ------------------------------------------------------------------------------
## .Qmin2 =
## function(etatry)
## { # A function implemented by Diethelm Wuertz
## # Description:
## # Internal Function
## # FUNCTION:
## # Compute:
## ans = .Qeta(etatry)$B
## assign("bBb", ans, pos = 1)
## # Return Value:
## ans
## }
## # ------------------------------------------------------------------------------
## # Internal Function: WHITTLE ESTIMATOR:
## .whittle =
## function(xinput, nsub = 1, model = c("farma", "fgn"), pp = 1,
## qq = 1, h = 0.5, ar = c(0.5), ma = c(0.5), out = TRUE, spec = FALSE)
## { # A function implemented by Diethelm Wuertz
## # Description:
## # Internal Function
## # Author:
## # Jan Beran; modified: Martin Maechler, Diethelm Wuertz 2004
## # FUNCTION:
## # To suppress No Visible Bindings Warning/Error:
## if(FALSE) { n = p = q = imodel = out = yper = bBb = NA }
## # Settings:
## model = model[1]
## assign("out", out, pos = 1)
## nmax = length(xinput)
## startend = c(1, nmax)
## istart = startend[1]
## iend = startend[2]
## nloop = nsub
## assign("n", trunc((iend - istart+1)/nloop), pos = 1)
## nhalfm = trunc((n - 1)/2)
## if(model == "farma")
## assign("imodel", 2, pos = 1)
## if(model == "fgn")
## assign("imodel", 1, pos = 1)
## assign("p", 0, pos = 1)
## assign("q", 0, pos = 1)
## p <- q <- 0
## if(imodel == 2) {
## assign("p", pp, pos = 1); assign("q", qq, pos = 1)
## } else {
## assign("p", 0, pos=1); assign("q", 0, pos = 1)
## }
## eta = c(h)
## if(p > 0) eta[2:(p+1)] = ar
## if(q > 0) eta[(p+2):(p+q+1)] = ma
## M = length(eta) #loop
## thetavector = c()
## i0 = istart
## flax = vector("list", nloop)
## # Necessary to make nsub/nloop work. VT.
## for (iloop in (1:nloop)) {
## h = max(0.2, min(h, 0.9))
## eta[1] = h
## i1 = i0+n - 1
## y = xinput[i0:i1]
## # Standardize Data:
## vary = var(y)
## y = (y - mean(y))/sqrt(var(y))
## # Periodogram of the Data:
## if(spec) {
## assign("yper", .per(y)[2:(nhalfm+1)], pos = 1)
## } else {
## assign("yper", .per(y)[2:(nhalfm+1)], pos = 1)
## }
## s = 2*(1-h)
## etatry = eta
## # Modified to make optim not give incorrect result. VT
## if(imodel == 1) {
## result = optim(par = etatry, fn = .Qmin2,
## method = "L-BFGS-B", lower = 0, upper = 0.999)
## } else {
## result = optim(par = etatry, fn = .Qmin2)
## }
## eta = result$par
## sturno = result$message
## theta1 = .Qeta(eta)$theta1
## theta = c(theta1, eta)
## thetavector = c(thetavector, theta)
## # Calculate goodness of fit statistic
## Qresult = .Qeta(eta) #output
## M = length(eta)
## if(imodel == 1) {
## SD = .CetaFGN(eta)
## SD = matrix(SD, ncol = M, nrow = M, byrow = TRUE) / n
## } else {
## # Changed to eliminate crashing in solve.qr in CetaARIMA. VT
## cat("M =", M, "\n")
## if(M > 2) {
## for (i in 3:M) {
## for (j in 2:(i-1)) {
## temp = eta[i]+eta[j]
## if(abs(temp) < 0.0001) {
## cat("Problem with estimating confidence intervals,",
## "parameter ", i, "and parameter ", j,
## "are the same, eliminating.\n")
## eta = eta[ - i]
## eta = eta[ - j]
## M = M - 2
## p = p - 1
## q = q - 1
## }
## }
## }
## }
## SD = .CetaARIMA(eta, p, q)
## SD = matrix(SD, ncol = M, nrow = M, byrow = TRUE)/n
## }
## Hlow = eta[1] - 1.96 * sqrt(SD[1, 1])
## Hup = eta[1]+1.96 * sqrt(SD[1, 1])
## if(out) {
## cat("theta =", theta, fill = TRUE)
## cat("H =", eta[1], fill = TRUE)
## cat("95%-CI for H: [", Hlow, ",", Hup, "]", fill = TRUE)
## }
## # Changing of the signs of the moving average parameters
## # in order to respect the sign of the Splus convention
## if(q > 0) eta[(p+2):(p+q+1)] = -eta[(p+2):(p+q+1)]
## etalow = c()
## etaup = c()
## for (i in (1:length(eta))) {
## etalow = c(etalow, eta[i] - 1.96 * sqrt(SD[i, i]))
## etaup = c(etaup, eta[i]+1.96 * sqrt(SD[i, i]))
## }
## if(out) {
## cat("95%-CI:", fill = TRUE)
## print(cbind(etalow, etaup), fill = TRUE)
## }
## if(spec) {
## cat("Periodogram is in yper", fill = TRUE)
## assign("fest", Qresult$theta1 * Qresult$fspec, pos = 1)
## cat("Spectral density is in fest", fill = TRUE)
## }
## flax[[iloop]] = list()
## flax[[iloop]]$par = eta
## flax[[iloop]]$sigma2 = bBb * var(xinput)
## flax[[iloop]]$conv.type = sturno
## remove("bBb", pos = 1)
## # Next subseries:
## i0 = i0+n
## # Changing of the signs of the moving average parameters:
## if(q > 0) eta[(p+2):(p+q+1)] = -eta[(p+2):(p+q+1)]
## } # end of nloop
## # Return:
## return(flax)
## # Clean up:
## remove("n", pos = 1)
## remove("p", pos = 1)
## remove("q", pos = 1)
## remove("imodel", pos = 1)
## remove("out", pos = 1)
## if(spec == FALSE) remove("yper", pos = 1)
## }
## ################################################################################
Try the fArma package in your browser
Any scripts or data that you put into this service are public.
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
Embedding an R snippet on your website
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.