library(bspec)
library(lomb)
library(WaveletComp)
library(zoo)
library(quantmod)
library(dplyr)
library(astrochron)
library(IRISSeismic)
library(Hmisc)
setwd("D:/Masters in Climate Change/Spectral Analysis/Data")

1.1

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

1.4

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

1.5

1/0.095 gives 10.52632 which corresponds to the cycles within the data.

1.6

rm(spec1)

2.1

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 Carrying out a Welch spectral analysis.

Naowelch <- welchPSD(NAO,seglength = 100)
plot(Naowelch$frequency,Naowelch$power,type="l",xlim=c(0,0.6))

2.3

Naowelch$frequency[findPeaks(Naowelch$power, thresh=1)] 
## [1] 0.09 0.14 0.21
These are the 3 largest peaks 0.09 0.14 0.21

2.4

1/0.09 gives 11.11111, 1/0.14 gives 7.142857, 1/0.21 gives 4.761905,

2.5

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

3.1

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

3.2

Ireprecip<- ts(mydata,start=1825, end=2017)
plot (Ireprecip, type="l", main='Irish precipitation timeseries')

3.3

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

3.4

data_winter <- aggregate(data2[,3],list(data2$Year), mean) 
detrended <- detrend(data_winter) 
## 
## ----- SUBTRACTING LINEAR TREND FROM STRATIGRAPHIC SERIES -----
##  * Slope= 0.07684393 
##  * y-intercept= -51.38379

dev.off()
## null device 
##           1
detrended_ts <- ts(detrended$x,start=1711, end=2017)
plot (detrended_ts, type="l", main='Detrended Irish precipitation timeseries')

3.6

detrended_tswelch <- welchPSD(detrended_ts,seglength = 50)
plot(detrended_tswelch$frequency,detrended_tswelch$power,type="l",xlim=c(0,0.6))

3.7

detrended_tswelch$frequency[findPeaks(detrended_tswelch$power, thresh=600)] 
## [1] 0.04 0.14 0.20 0.30
the four largest peaks were 0.04 0.14 0.20 0.30

3.8

1/0.04 gives 25
1/0.14 gives 7.142857
1/0.20 gives 5
1/0.30 gives 3.333333

3.9

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

4.1

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

ts1<- ts(precipitation_winter$V1,start=1711, end=2016) 
ts2<- ts(NAO_winter$V2,start=1825, end=2017)
this turns it into a timeseries record

4.3

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

4.4

correlation = rcorr(t$ts1,t$ts2, type="pearson")
The two time series correlate significantly based on the correlation table. The result of 0.12 significance is larger then 0.01 significance tested againts.

4.5

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

4.6

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

4.7

4.8

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

5.1

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

5.2

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

### There are very little cycles to examine, however there are peaks around each interval of 1000 followed by a decrease. ### The amplitude of each peak differs however there is arguably a cylce in the Age axis.

5.3

mydata$Depth=NULL

5.4

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

5.5

The signinficant cycles are at 0.001 and 0.003 frequency and there are only two peaks.

5.6

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

6.1

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

6.2

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

6.3

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)

6.4

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

7.1

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

7.2

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

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

There is one main signiicant cycle that occurs between 2500-3500 years ago. The significant level is approximately 500.

There are other cycles over periods much smaller and they occur more freuently, however they are less significant.

7.5

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

8.1

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

8.2

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

length is the length of NAO_precip_winter

8.3

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

8.4

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

8.5

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

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

Note that the echo = FALSE parameter was added to the code chunk to prevent printing of the R code that generated the plot.