\documentclass[article,nojss]{jss} %% NOTA BENE: More definitions --> further down %%%%%%%%%%%% % \author{Martin M\"achler \\ ETH Zurich% \\ June 2011 {\tiny (\LaTeX'ed \today)}%---- for now } \title{Numerical Explorations for Fast Spectrum of Fractional Gaussian Noise} % \def\mythanks{a version of this paper, for \pkg{nacopula} 0.4\_4, has been published % in JSS, \url{http://www.jstatsoft.org/v39/i09}.} %% for pretty printing and a nice hypersummary also set: \Plainauthor{Martin M\"achler} %% comma-separated %\Plaintitle{Numerical Explorations for Fast Spectrum of Fractional Gaussian Noise} \Shorttitle{Fast Spectrum of Fractional Gaussian Noise} % %\VignetteIndexEntry{fast-specFGN} %\VignetteDepends{longmemo} %\VignetteDepends{sfsmisc} \SweaveOpts{engine=R,eps=FALSE,pdf=TRUE,width=7,height=5,strip.white=true,keep.source=TRUE} %% an abstract and keywords \Abstract{ The package \pkg{longmemo} .... Paxson (1997) ... ... } \Keywords{Euler-Maclaurin Formula, Fractional Gaussian Noise, Spectrum} %% at least one keyword must be supplied %% publication information %% NOTE: Typically, this can be left commented and will be filled out by the technical editor %% \Volume{13} %% \Issue{9} %% \Month{September} %% \Year{2004} %% \Submitdate{2004-09-29} %% \Acceptdate{2004-09-29} %% The address of (at least) one author should be given %% in the following format: \Address{ Martin M\"achler\\ Seminar f\"ur Statistik, HG G~16\\ ETH Zurich\\ 8092 Zurich, Switzerland\\ E-mail: \email{maechler@stat.math.ethz.ch}\\ URL: \url{http://stat.ethz.ch/people/maechler} } %% It is also possible to add a telephone and fax number %% before the e-mail in the following format: %% Telephone: +43/1/31336-5053 %% Fax: +43/1/31336-734 %% for those who use Sweave please include the following line (with % symbols): %% MM: this is "substituted" by jss.cls: %% need no \usepackage{Sweave.sty} \usepackage[american]{babel}%for American English \usepackage{amsmath}%sophisticated mathematical formulas with amstex (includes \text{}) \usepackage{mathtools}%fix amsmath deficiencies \usepackage{amssymb}%sophisticated mathematical symbols with amstex (includes \mathbb{}) % \usepackage{amsthm}%theorem environments \usepackage{bm}%for bold math symbols: \bm (= bold math) \usepackage{enumitem}%for automatic numbering of new enumerate environments % This is already in jss above -- but withOUT the fontsize=\small part !! \DefineVerbatimEnvironment{Sinput}{Verbatim}{fontsize=\small,fontshape=sl} \DefineVerbatimEnvironment{Soutput}{Verbatim}{fontsize=\small} \DefineVerbatimEnvironment{Scode}{Verbatim}{fontsize=\small,fontshape=sl} %% makes space between Sinput and Soutput smaller: \fvset{listparameters={\setlength{\topsep}{0pt}}}% !! quite an effect! %% ??? FIXME but it also influences all other lists, itemize, ... ???? FIXME %% \setkeys{Gin}{width=\textwidth}% Sweave.sty has {width=0.8\textwidth} \newcommand*{\R}{\proglang{R}}%{\textsf{R}} \newcommand*{\eps}{\varepsilon} \newcommand{\abs}[1]{\left| #1 \right|}% \abs{ab} --> | ab | ``absolut Betrag'' \newcommand{\norm}[1]{\left\| #1 \right\|}%- \norm{ab} --> || ab || \newcommand*{\BB}{\mathcal{B}(\lambda, H)} %% journal specific aliases \newcommand*{\setcapwidth}[1]{} %% end of declarations %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% \begin{document} %% include your article here, just as usual %% Note that you should use the \pkg{}, \proglang{} and \code{} commands. % \section[About Java]{About \proglang{Java}} %% Note: If there is markup in \(sub)section, then it has to be escape as above. %% Note: These are explained in '?RweaveLatex' : <>= op.orig <- options(width = 75, SweaveHooks= list(fig=function() par(mar=c(5.1, 4.1, 1.1, 2.1))), useFancyQuotes = FALSE, ## for JSS, but otherwise MM does not like it: ## prompt="R> ", continue=" ")# 2 (or 3) blanks: use same length as 'prompt' Sys.setenv(LANGUAGE = "en") if(.Platform$OS.type != "windows") Sys.setlocale("LC_MESSAGES","C") @ %\section[Introduction]{Introduction \small~\footnote{\mythanks}} %\section{Introduction} \section{.. intro ..} The spectral density of fractional Gaussian noise (``fGn'') with Hurst parameter $H \in (0,1)$ is (\cite{BerJ86, BerJ94}) \begin{equation} \label{spec-def} f_H(\lambda) = \mathcal{A}(\lambda, H) \left(\abs{\lambda}^{-2H-1} + \BB \right), \end{equation} for $\lambda\in [-\pi, \pi]$, where $\mathcal{A}(\lambda, H) == 2\sin(\pi H)\, \Gamma(2H+1)\, (1 - \cos\lambda)$, and \begin{equation} \label{def-B} \BB = \sum_{j=1}^\infty \left( (2\pi j + \lambda)^{-(2H+1)} + (2\pi j - \lambda)^{-(2H+1)} + \right). \end{equation} For the Whittle estimator of $H$ and also other purposes, its advantageous to be able to evaluate $f_H(\lambda_i)$ efficiently for a whole vector of $\lambda_i$, typically Fourier frequencies $\lambda_i = 2\pi i/n$, for $i=1,2,\dots,\lfloor (n-1)/2 \rfloor$. Such evaluation is problematic because of the infinite sum for $\BB$ in (\ref{def-B}). Traditionally, e.g., already in Appendix.... of \citet{BerJ94}, the infinite sum $\sum_{j=1}^\infty$ had been replaced by $\sum_j^{200}$ --- which was still not very efficient and not extremely accurate. In our \R\ package \pkg{longmemo}, we now provide the function \code{B.specFGN}$(\lambda$, $H)$ to compute $\BB$, using several ways to compute the infinite sum approximately, e.g., for $H = 0.75$ and $n = 500$, i.e., at 250 Fourier frequencies, <>= require("longmemo") fr <- .ffreq(500) B.1 <- B.specFGN(fr, H = 0.75, nsum = 200, k.approx=NA) B.xct <- B.specFGN(fr, H = 0.75, nsum = 10000, k.approx=NA) all.equal(B.xct, B.1) @ which means that the 200 term approximation is accurate to \ 4 \ decimal digits for $H = .75$ but the accuracy is smaller for smaller $H$. %% FIXME: give numbers %% use B.specFGN() For this reason, \citet{PaxsonV:1997:FGN} derived formulas for fast and stilly quite accurate approximations of $\BB$, noting that $\BB = \sum_{j=1}^\infty f(j; \lambda,H)$ for \begin{equation} \label{def-f} f(x; \lambda,H) = (2\pi x + \lambda)^{-(2H+1)} + (2\pi x - \lambda)^{-(2H+1)}, \end{equation} and the fact that $\sum_{j=1}^\infty f(j)$ is a Riemann sum approximation of $\int_0^\infty f(x)\;\mathrm{d}x$ or $\int_1^\infty f(x)\;\mathrm{d}x$. <>= ## ^^^^^^ drop comments here ##' @title The function f(x) used in the approximation of B.specFGN() ##' @param x numeric vector of positive values ##' @param lambda number in [-pi, pi]; practically in (0, pi] ##' @param H Hurst parameter in (0, 1); practically in (1/2, 1) ##' @return ##' @author Martin Maechler fB <- function(x, lambda, H) { u <- 2*pi*x h <- -(2*H+1) (u + lambda)^h + (u - lambda)^h } @ Now its clear that $f(x)$ cannot be computed (or ``is infinite'') at $x=0$, and more specifically, $f(x)$ tends to $\infty$ when $x\to \frac{\lambda}{2\pi}$, as in the second term of $f$, $2\pi x - \lambda$ only remains positive when $2\pi x > \lambda$. This is always fulfilled for $x in \{1,2,\dots\}$, as $\lambda < \pi$, but is problematic when considering $\int_0^b f(x)\;\mathrm{d}x$ as above. % Some illustrations of the function $f(x;\lambda, H)$ and its ``pole'' at $\frac{\lambda}{2\pi}$ : <>= draw.f <- function(lambda, H, ..., log = "", ylim= if(any("y" == strsplit(log,"")[[1]]))NULL else c(0, 2*fB1)) { fB1 <- fB(1, lambda=lambda, H = H) curve(fB(x, lambda=lambda, H = H), ylim=ylim, log=log, xaxt="n", ...) sfsmisc::eaxis(1) mtext(bquote(list(lambda == .(lambda), H == .(H)))) abline(v = lambda / (2*pi), lty=2, lwd=3, col = "blue3") xrng <- par("usr")[1:2]; if(xlog <- par("xlog")) xrng <- 10^xrng abline(h = 0, lty=3, col="gray20") j <- floor(xrng[1]):ceiling(xrng[2]) lines(j, fB(j, lambda=lambda, H = H), type = "h", lty=5, lwd= 0.75, col = "gray40") axis(1, at=j[j >= if(xlog) 2 else 1][1:2], col="green3") } op <- sfsmisc::mult.fig(4)$old.par draw.f(lambda=.001, H = 0.75, 0, 10) draw.f(lambda= 0.5, H = 0.75, 0, 10) draw.f(lambda= 1.5, H = 0.75, 0, 10) draw.f(lambda= 3.0, H = 0.75, 0, 10) <>= sfsmisc::mult.fig(3)$old.par draw.f(lambda=.001, H = 0.99, .001, 100, log="xy") draw.f(lambda= .5, H = 0.99, .001, 100, log="xy") ## ! draw.f(lambda= 3.0, H = 0.51, .001, 100, log="xy") ## ! par(op) @ So, very clearly, Paxson's first formula, using $\int_0^1 f(x)\; \mathrm{d}x$ is not feasible, as $f(x)$ is \emph{not} defined (or defined as $\infty$) for $x \le \lambda/(2\pi)$. However, his generalized formula, ``(7), p.~15'', \begin{equation} \label{eq:Paxson-2} \sum_{i=1}^\infty f_i \approx \sum_{j=1}^k f_j + \frac 1 2\int_k^{k+1} f(x)\;\mathrm{d}x + \int_{k+1}^\infty f(x)\;\mathrm{d}x, \end{equation} clearly \emph{is} usable for $k \ge 1$ (but not for $k=0$, contrary to what he suggests). Indeed, with \code{B.specFGN(}{$\lambda, H$,} \code{k.approx)}, we now provide the result of applying approximation (\ref{eq:Paxson-2}) to the infinite sum for $\BB$ in (\ref{def-B}). Paxson ended the $k=3$ approximation which he further improved considerably, empirically, by numerical comparison (and least squares fitting) with the ``accurate'' formula using \code{nsum = }$10'000$ terms. In the following section, we propose another improvement over Paxson's original idea: \section{Better approximations using the Euler--Maclaurin formula} Copied straight from %% FIXME !! {\large \url{http://en.wikipedia.org/wiki/Euler-Maclaurin_formula}} : If $n$ is a natural number and $f(x)$ is a smooth, i.e., sufficiently often differentiable function defined for all real numbers $x$ between 0 and $n$, then the integral \begin{equation} \label{eq:I} I=\int_0^n f(x)\,dx \end{equation} can be approximated by the sum (or vice versa) \[ S=\frac{1}{2}f(0)+ f(1) +\cdots+ f(n-1) +\frac{1}{2}f(n) \] (see trapezoidal rule). The Euler--Maclaurin formula provides expressions for the difference between the sum and the integral in terms of the higher derivatives $f^{(k)}$ at the end points of the interval 0 and n. Explicitly, for any natural number $p$, we have \[ S-I = \sum_{k=2}^p {\frac{B_k}{k!}\left(f^{(k-1)}(n) - f^{(k-1)}(0)\right)} + R \] where $B_1 = -1/2$, $B_2 = 1/6$, $B_3 = 0$, $B_4 = -1/30$, $B_6 = 1/42$, $B_8 = -1/30$, \dots are the Bernoulli numbers, and $R$ is an error term which is normally small for suitable values of p. (The formula is often written with the subscript taking only even values, since the odd Bernoulli numbers are zero except for $B_1$.) Note that \[ -B_1(f(n)+f(0)) =\frac{1}{2}(f(n)+f(0)). \] Hence, we may also write the formula as follows: \begin{equation} \label{eq:EM-form} \sum_{i=0}^n f(i) = \int^n_0f(x)\,dx-B_1(f(n)+f(0))+\sum_{k=1}^p\frac{B_{2k}}{(2k)!}\left(f^{(2k-1)}(n)-f^{(2k-1)}(0)\right)+R. \end{equation} ............. In the context of computing asymptotic expansions of sums and series, usually the most useful form of the Euler--Maclaurin formula is \[ \sum_{n=a}^b f(n) \sim \int_a^b f(x)\,dx + \frac{f(a)+f(b)}{2} + \sum_{k=1}^\infty \,\frac{B_{2k}}{(2k)!}\left(f^{(2k-1)}(b)-f^{(2k-1)}(a)\right), \, \] where $a$ and $b$ are integers. %[2] Often the expansion remains valid even after taking the limits ${a\to -\infty}$ or ${b\to +\infty}$, or both. % In many cases the integral on the right-hand side can be evaluated in closed form in terms of elementary functions even though the sum on the left-hand side cannot. {\footnotesize (end of citation from Wikipedia)} --- --- \section{Session Information} <>= toLatex(sessionInfo()) <>= options(op.orig) @ \section{Conclusion} \bibliography{longmemo} \end{document}