library(bspec)
library(lomb)
library(WaveletComp)
library(zoo)
library(quantmod)
library(dplyr)
library(astrochron)
library(IRISSeismic)
library(Hmisc)
setwd("C:/MSc_CC/lisas")

1.1: ‘Sunspots’ is a dataset available in R already - plot this data

1.2: Compute and plot the Welch method spectrum:

spec1 <- welchPSD(sunspots, seglength=200) 
plot(spec1$frequency, spec1$power, log='y',type="l", xlim=c(0,1))

1.4: What are the frequencies of the peaks?

spec1$frequency[findPeaks(spec1$power, thresh=30000.0)] 
## [1] 0.095

1.5: What are the cycles within the data?

1/0.095
## [1] 10.52632

1.6 Removing unwanted files

rm(spec1)

North Atlantic Oscillation data

2.1: Import the NOA data and view the record

read.csv("NAO.csv", header=TRUE)
mydata <- read.table("NAO.csv", header=FALSE, sep=",")[,2] 
NAO <- ts(mydata,start=1825, end=2017)
plot (NAO, type="l")

2.2 Carry out the welch analysis

NAOwelch <- welchPSD(NAO, seglength=100) # use seg100 for NAO
plot(NAOwelch$frequency, NAOwelch$power,type="l", xlim=c(0,0.6))

2.3 What are the cycles shown by the three largest peaks

NAOwelch$frequency[findPeaks(NAOwelch$power, thresh=1)]
## [1] 0.09 0.14 0.21

Peak 1

1/0.09
## [1] 11.11111

Peak 2

1/0.14
## [1] 7.142857

Peak 3

1/0.21
## [1] 4.761905

2.4 Export the data

write.table(NAO, file = "Winter_NAO_data.csv", col.names=FALSE, sep=",")

2.5 Clear environment

rm(NAOwelch)
rm(mydata)
rm(NAO)

3.1 Import the data ‘Precipitation_Ireland.csv’

read.csv("Precipitation_Ireland.csv", header=TRUE) 
mydata <- read.table("Precipitation_Ireland.csv", header=TRUE, sep=",")

3.2 Plot the data

IrePrec <- ts(mydata,start=1711, end=2016) 
plot(IrePrec, type="l")

3.3 Select the rainfall from Jan and Feb

data2 <- (mydata %>% filter(mydata$Month == c(1, 2)))

3.4 Extract average rainfall for each year for Jan and Feb

data_winter <- aggregate(data2[,3],list(data2$Year), mean)

3.5 Detrend

detrended <- detrend(data_winter) 

dev.off()

3.6 Convert to a timeseries data and plot

detrendedts <- ts(detrended$x,start=1711, end=2016) 
plot(detrendedts, type= "l")

3.6 Perform the Welch method

detrendtswelch <- welchPSD(detrendedts, seglength=50) # use seg50 for ire prec
plot(detrendtswelch$frequency, detrendtswelch$power,type="l", xlim=c(0,0.6))

3.7 what are the cycles- threshold 600

detrendtswelch$frequency[findPeaks(detrendtswelch$power, thresh=600)]
## [1] 0.04 0.14 0.20 0.30

Peak 1

1/0.04
## [1] 25

Peak 2

1/0.14
## [1] 7.142857

Peak 3

1/0.2
## [1] 5

Peak 4

1/0.3
## [1] 3.333333

3.8 Export

write.table(detrendedts, file = "Winter _precipitation_data.csv", row.names=FALSE, col.names=FALSE, sep=",")

Part B: Cross-spectral Analysis

4.1 Install the detrended data on the NAO and winter precipitation from Part A

read.csv("Winter _precipitation_data.csv", header=FALSE) 
precipitation_winter <- read.table("Winter _precipitation_data.csv", sep=",")
read.csv("Winter_NAO_data.csv", header=FALSE) 
NAO_winter <- read.table("winter_NAO_data.csv", sep=",")

4.2 Convert each dataset to timeseries data.

ts1 <- ts(precipitation_winter$V1,start=1711, end=2016)
ts2 <- ts(NAO_winter$V2, start=1825, end=2017)

4.3 Combine into one dataframe

t <- ts.intersect(ts1,ts2,dframe=TRUE)

4.4 Using the table below assess whether the two timeseries correlate significantly?

correlation = rcorr(t$ts1,t$ts2, type="pearson")

The correlation value between winter precipitation and winter NAO is 0.12. The number of data points across both time series is 499. Reading across the Pearson correlation table we arrive at a value of 0.01 along the top. The correlation value is 0.12 is larger than this and therefore can be considered significant.

4.5 Calculate the cross spectral analysis between the two timeseries

DF <- crossSpectrum(t) 
transferFunction <- DF$Pxy / DF$Pxx 
transferAmp <- Mod(transferFunction)

4.6 Plot the results

plot(DF$freq,transferAmp,type='l',log='x',      xlab="Frequency",      main="Transfer Function Amplitude")

4.7 Which cycle is a feature of both datasets?

A clear cycle for both data sets can be seen at a frequency of 0.2. This reaches a peak transfer amplitude of 1.

4.8 Export the data, this will be used later

write.table(t, file = "NAO_and_precipitation_data.csv", col.names=FALSE, sep=",")

Part C: Unevenly spaced timeseries data

5.1 Load the data

read.csv("pedrido.csv", header=TRUE) 
mydata <- read.table("pedrido.csv", header=TRUE, sep=",")

5.2 Plot the data - can you see any cycles with a visual assessment?

plot(mydata$Sand~mydata$Age, type= "l")

There is possibly a cycle occurring at every 1000 intervals of age. Roughly, every interval of 1000 along the age axis there is a sudden increase to a peak in the sand data. This tends to be followed by a sudden decrease in the sand data. These peaks increase to different amplitudes, but there is roughly a 1000 age interval year frequency.

5.3 Remove column about depth

mydata$Depth=NULL

5.4 Using the “lomb” package

ls<- lsp(mydata, from = NULL, to = NULL, ofac = 5, alpha = 0.05, plot = TRUE)

5.5 What are the significant cycles?

A significance of over 95% (0.05) was chosen. There are two peaks at a frequency of 0.001 and 0.003.

5.6 Another approach is to plot as Power vs Period (rather than frequency)

lsp(mydata, from = NULL, to = NULL, ofac = 5, alpha = 0.05,type = c("period"), plot = TRUE)

6.1 Define the new age axis which is evenly spaced

New_age <- seq(from = -10, to = 4600, by = 10)

6.2 Interpolate and resample the data evenly - linear interpolation has been used

new_data <- approx(x=mydata$Age, y=mydata$Sand, xout=New_age, method="linear") 
names(new_data) <- c("date","sand")

6.3 Plot the original data against the downsampled data, how do they differ?

plot(mydata$Age, mydata$Sand, type="o", col="blue", lty=1, main="Original Data vs Downsampled Data") 
points(new_data$date, new_data$sand, col="red", pch="*") 
lines(new_data$date, new_data$sand, col="red",lty=2)

The original data appear to have occasionally more dramatic and higher peaks than the downsampled data. For example, at the age interval of around 900, the original sand data reaches a peak value of 12, while the downsampled data only reaches a peak of just over 10. Similarly, the peak at the age interval of 3300 shows a difference in peak values, with the original data peak reaching 15, and the downsampled data reaching 14. Generally, the original data and downscaled data follow the same pattern, with an occasional reduction in peak values in the downsampled data.

6.4 Export the evenly spaced data

write.csv(new_data, file = "PedridoEvenData.csv")

Part D

7.1 Importing and formatting data ready for using in wavelet transform

Pedrido <- read.csv("PedridoEvenData.csv")[ ,3] 
New_age <- read.csv("PedridoEvenData.csv")[ ,2] 
my.data = data.frame(date = New_age, Pedrido=Pedrido)

7.2 Wavelet transform analysis

my.wt = analyze.wavelet(my.data, "Pedrido", loess.span = 0, dt = 10, dj = 1/20, make.pval = TRUE, method = "white.noise", params = NULL, n.sim = 100, date.format = "%Y", date.tz = NULL, verbose = TRUE)

7.3 Plotting the wavelet transform

wt.image(my.wt, color.key = "i", label.time.axis = TRUE, main = "wavelet power spectrum", legend.params = list(lab = "wavelet power levels"), periodlab = "period (years)", timelab = "time (years)", date.format = "%d", spec.time.axis = list(at = seq(-10, 462, by = 10), labels = seq(0, 460, by = 50)))
## Warning in wt.image(my.wt, color.key = "i", label.time.axis = TRUE, main = "wavelet power spectrum", : 
## Please check your time axis specifications. Default settings were used.

7.4 Which cycles are significant? Over which time period do the cycles occur?

There appears to be one main significant cycle. This occurs between 2000-3500 year ago. This had a cyclical frequency of roughly 500 years. Two other cycles appear but seem to be relatively less significant.

7.5 Calculate the average cycles

wt.avg(my.wt, siglvl=0.05, sigcol="red")

Part E

8.1 Import the data of the NAO and precipitation

read.csv("NAO_and_precipitation_data.csv", header=FALSE) 
NAO_precip_winter <- read.table("NAO_and_precipitation_data.csv", sep=",")[,2:3]

8.2 Define the calender ages

my.date <- seq(as.POSIXct("1825-01-30 00:00:00", format = "%Y"), by = "year", length.out = 192)

8.3 Combine the new age column (my.date) with the NAO and precipitation data

date<- my.date 
my.data <- cbind(NAO_precip_winter,date)

8.4 Conduct the wavelet coherency analysis

my.wc <- analyze.coherency(my.data, c(1, 2), dt=1)

8.5 Plot the wavelet power spectrum

wc.image(my.wc, main = "cross-wavelet power spectrum, x over y", legend.params = list(lab = "cross-wavelet power levels"),          periodlab = "period (years)", show.date = TRUE)

8.6 Calculate the average cycles. What are the significant cycles shared by the precipitation and NAO data? (significant cycles are those marked by red points)

wc.avg(my.wc, siglvl = 0.05, sigcol = 'red', periodlab = "period (years)")

A significant cyclical event occurs every 8 years for both Irish Precip and NOA.

Part F

Reading in the Central England Temperature data

read.csv("CET.csv", header=TRUE)

Importing and formatting data ready for using in wavelet transform

CET <- read.csv("CET.csv")[ ,14] 
Year <- read.csv("CET.csv")[ ,1] 
my.data = data.frame(date = Year, CET=CET)
Wavelet transform analysis
my.wt = analyze.wavelet(my.data, "CET", loess.span = 0, dt = 10, dj = 1/20, make.pval = TRUE, method = "white.noise", params = NULL, n.sim = 100, date.format = "%Y", date.tz = NULL, verbose = TRUE)

Plotting the wavelet transform

wt.image(my.wt, color.key = "i", label.time.axis = TRUE, main = "wavelet power spectrum", legend.params = list(lab = "wavelet power levels"), periodlab = "period (years)", timelab = "time (years)", date.format = "%d", spec.time.axis = list(at = seq(-10, 462, by = 10), labels = seq(0, 460, by = 50)))
## Warning in wt.image(my.wt, color.key = "i", label.time.axis = TRUE, main = "wavelet power spectrum", : 
## Please check your time axis specifications. Default settings were used.

Which cycles are significant? Over which time period do the cycles occur?

There seems to be less pronounced significant events in this power spectrum. The most significant events occurs around 2750 years ago. This had a cyclical frequency of around 16 years. Another significant event seems to have occurred around around 2300-2400 years ago. This had a cyclical frequency of around 16 years also.

Calculate the average cycles

wt.avg(my.wt, siglvl=0.05, sigcol="red")