simon_example.R

# you will need my package "gtseries"
# for this you  first need to install the 'devtools' standard R package, and then: 

# the lines below should do this for you. If it does not work, install the "devtool" standard package, 
# and then type
# devtools::install_git('https://forge.uclouvain.be/mcrucifix/tseries.git')


list.of.gitpackages <- c("gtseries")
new.gitpackages <- list.of.gitpackages[!(list.of.gitpackages %in% installed.packages()[,"Package"])]
if (length(new.gitpackages)) { 
  list.of.packages <- c("devtools")
  few.packages <- list.of.packages[!(list.of.packages %in% installed.packages()[,"Package"])]
  if (length(new.packages)) install.packages(new.packages)
  devtools::install_git('https://forge.uclouvain.be/mcrucifix/tseries.git')
}




data <- read.csv('example Aalst agri.csv')
temp <- data$temp[-1]

# note: the below does not dela with Na
temp <- ts(temp,deltat=0.5)

require(gtseries)

C = cwt_morlet(temp)
C
plot(C)

range.fast <- which(attr(C,"scale") < 20.0 )
range.slow <- which(attr(C,"scale") > 36 )
range.daily  <- seq(length(attr(C,"scale")))[-c(range.fast,range.slow)]

signal.fast = reconstruct_morlet(C,scales=c(0,12))
signal.slow = reconstruct_morlet(C,scales=c(60,Inf))
signal.daily = reconstruct_morlet(C,scales=c(12,60))


signal.total  <-  Re(signal.slow) + Re(signal.fast) + Re(signal.daily)

plot(temp - mean(temp), type='l')
lines(ts(Re(signal.fast), deltat=.5), col='red', type='l')  # fast variability, in red
lines(Re(signal.daily), col='blue', type='l')               # daily variability, in blue
lines(Re(signal.slow), col='green', type='l')               # slow part (trend), in green
lines(Re(signal.total), col='purple', type='l')             # reconstruction (not perfecct)

# The Amplitude of the daily cycle can be obtained as follows, but it is surprisingly smooth and high-frequency, so probably it is best diagnosed by hand.

Amp.daily  <-  Mod(signal.daily)

plot(signal.daily)                    # the daily signal
lines(Im(signal.daily), col='red')    # 90-degree phase lagged (in red)
lines(Mod(signal.daily), col='blue')  # amplitude, in Blue

# Similary, the amplitude of the fast variability may be obtained with: 

plot(Mod(signal.fast))                    # the daily signal

# though again the result in not quite convincing (much high frequency).