From 5247fe5ca083a4bba61bdd604e9efc7b3cd15913 Mon Sep 17 00:00:00 2001
From: Michel Crucifix <michel.crucifix@uclouvain.be>
Date: Wed, 22 May 2024 22:05:31 +0200
Subject: [PATCH] haar

---
 .Rbuildignore                                |  1 +
 R/haar.R                                     | 48 ++++++++++++++++++++
 miscellaneous_notes/fluctuation_spectrum.txt | 18 ++++++++
 src/Makevars                                 |  4 +-
 4 files changed, 69 insertions(+), 2 deletions(-)
 create mode 100644 R/haar.R
 create mode 100644 miscellaneous_notes/fluctuation_spectrum.txt

diff --git a/.Rbuildignore b/.Rbuildignore
index 6bbefad..40c7223 100644
--- a/.Rbuildignore
+++ b/.Rbuildignore
@@ -1,6 +1,7 @@
 .gitignore
 test
 old_man
+miscellaneous_notes
 ^data-raw$
 .Rprofile
 TODO
diff --git a/R/haar.R b/R/haar.R
new file mode 100644
index 0000000..52dcba9
--- /dev/null
+++ b/R/haar.R
@@ -0,0 +1,48 @@
+# TO DO: THERE IS A BIAS FOR THE LONG LENGHTSCALES
+# I GUESS MUST MULTIPLY BY (N-1)/N where N is 
+# length / scale
+
+
+
+#' Haar fluctuation spectrum
+#'
+#' @description Haar fluctuation spectrum after Lovejoy
+#' still needs to be documented
+#' currently requires n=2^i data points
+
+
+#' @param x the input series (typically a numeric)
+#' @param q the fluctuation order
+#' @examples
+#' x = rnorm(2048)
+#' xi1 = haar(x, q=1)
+#' xi2 = haar(x, q=2)
+
+haar <- function(x, q=2, discarded_scales = 4)
+{
+  n = length(x)
+  nsteps=floor(log(n)/log(2))
+  nmax = 2^nsteps
+  DT = sapply (seq(nsteps), function(i)
+  {
+     m = 2^i
+     nn = nmax / m
+     M = matrix(x, nmax/m, m, byrow=TRUE)
+     M1 = M[,1:(m/2), drop=FALSE]
+     M2 = M[,(m/2+1):m, drop=FALSE]
+     tM1 = apply(M1, 1, mean)
+     tM2 = apply(M2, 1, mean)
+     abs(tM2 - tM1) * ( 2 * nn ) / (2 * nn  - 1)
+  })
+  x= 2^seq(nsteps)
+  y= sapply(DT, function(x) mean(x^q)) 
+  H = data.frame(x = x,  y=y, logx=log(x), logy=log(y))
+  kept_scales <- seq(max(2,nsteps-discarded_scales))
+  out = lm( H$logy[kept_scales] ~ H$logx[kept_scales] )$coefficients[2]
+  attr(out, "scale") = H
+  attr(out, "q") = q
+  attr(out,"class") = "fluctuation spectrum"
+  out
+ }
+ 
+
diff --git a/miscellaneous_notes/fluctuation_spectrum.txt b/miscellaneous_notes/fluctuation_spectrum.txt
new file mode 100644
index 0000000..8ae1115
--- /dev/null
+++ b/miscellaneous_notes/fluctuation_spectrum.txt
@@ -0,0 +1,18 @@
+ligne rouge : cumsum de la simulation avec package 'FGN'
+ligne noire : simul avec FGN
+
+x-axis: le Hurst entre dans FGN
+y=axis: le H calcule par Haar Wavelet fluct. 
+
+a partir d'ici: on parle du H de Lovejoy
+
+H s'appelle H is the fluctua- tion (also called “nonconservation”) exponent 
+
+H=-0.5 = bruit blanc 
+
+H approx xi(1)
+H = xi(2) / 2 
+
+beta = 1 +  xi(2) = 1 + 2H
+
+donc, bruit blanc : xi(2) =  1 + 2*(-0.5) = 0
diff --git a/src/Makevars b/src/Makevars
index c2452b0..2871a83 100644
--- a/src/Makevars
+++ b/src/Makevars
@@ -1,7 +1,7 @@
 # Kindly supplied by Dirk Eddelbuettel
 # set by configure
-GSL_CFLAGS = -I/usr/local/include
-GSL_LIBS   = -L/usr/local/lib -lgsl -lgslcblas -lm
+GSL_CFLAGS = -I/usr/include
+GSL_LIBS   = -L/usr/lib/x86_64-linux-gnu -lgsl -lgslcblas -lm
 
 # combine to standard arguments for R
 PKG_CPPFLAGS =  $(GSL_CFLAGS) -I.
-- 
GitLab