library(ggplot2) library(gridExtra) www<-"https://www.mimuw.edu.pl/~noble/courses/TimeSeries/data/nino3data.asc" nino<-read.table(www,skip=3,header=F) names(nino) <- c("Year", "SST", "SSA") plot(nino$Year, nino$SST, type = "l") plot1 <- ggplot(data = nino) + geom_line(aes(y = SST, x = Year)) plot2 <- ggplot(data = nino) + geom_line(aes(y = SSA, x = Year)) grid.arrange(plot1, plot2) acf1 <- acf(nino$SST, lag.max = 12 * 20, plot = FALSE) acf2 <- acf(nino$SSA, lag.max = 12 * 20, plot = FALSE) plot1 <- ggplot() + geom_line(aes(x = acf1$lag/12, y = acf1$acf)) plot2 <- ggplot() + geom_line(aes(x = acf2$lag/12, y = acf2$acf)) grid.arrange(plot1, plot2) # Create dataframe with different harmonics X <- data.frame(Year=nino$Year, y = nino$SST, sin(2*pi*1*nino$Year), cos(2*pi*1* nino$Year), # sine and cos for frequency = 1 sin(2*pi*2*nino$Year), cos(2*pi*2*nino$Year), # freq. equals 2 (i.e. period= 6 months) sin(2*pi*1/3*nino$Year), cos(2*pi*1/3*nino$Year), # freq = 1/3 (period=3 years) sin(2*pi*1/3.5*nino$Year), cos(2*pi*1/3.5*nino$Year), # freq=3.5 (period=3.5 years) sin(2*pi*1/6*nino$Year), cos(2*pi*1/6*nino$Year), # freq=6 (period=6 years) sin(2*pi*1.01*nino$Year), cos(2*pi*1.01*nino$Year) # freq=1.01 (period=.99 years) ) ggplot(data=subset(X, Year>1980)) + geom_line(aes(x=Year, y=X[X$Year>1980,3])) ggplot(data=subset(X, Year>1980)) + geom_line(aes(x=Year, y=X[X$Year>1980,5])) ggplot(data=subset(X, Year>1980)) + geom_line(aes(x=Year, y=X[X$Year>1980,7])) ggplot(data=subset(X, Year>1980)) + geom_line(aes(x=Year, y=X[X$Year>1980,9])) ggplot(data=subset(X, Year>1980)) + geom_line(aes(x=Year, y=X[X$Year>1980,11])) mod <- lm(y ~ . - Year, data = X) # Regress y on everything (but Year) summary(mod) X$resid <- residuals(mod) X$pred <- predict(mod) ggplot(data = subset(X, Year > 1970)) + geom_line(aes(x = Year, y = y)) + geom_line(aes(x = Year, y = pred), color = "red") raw.spec <- spec.pgram(nino$SST, taper = 0) plot(raw.spec) plot(raw.spec, log = "no") # spec.df <- as.data.frame(raw.spec) spec.df <- data.frame(freq = raw.spec$freq, spec = raw.spec$spec) # Create a vector of periods to label on the graph, units are in years yrs.period <- rev(c(1/6, 1/5, 1/4, 1/3, 0.5, 1, 3, 5, 10, 100)) yrs.labels <- rev(c("1/6", "1/5", "1/4", "1/3", "1/2", "1", "3", "5", "10", "100")) yrs.freqs <- 1/yrs.period * 1/12 #Convert annual period to annual freq, and then to monthly freq spec.df$period <- 1/spec.df$freq ggplot(data = subset(spec.df)) + geom_line(aes(x = freq, y = spec)) + scale_x_continuous("Period (years)", breaks = yrs.freqs, labels = yrs.labels) + scale_y_continuous() ggplot(data = subset(spec.df)) + geom_line(aes(x = freq, y = spec)) + scale_x_continuous("Period (years)", breaks = yrs.freqs, labels = yrs.labels) + scale_y_log10() ggplot(data = subset(spec.df)) + geom_line(aes(x = freq, y = spec)) + scale_x_log10("Period (years)", breaks = yrs.freqs, labels = yrs.labels) + scale_y_log10() plot(kernel("daniell", m = 10)) # A short moving average plot(kernel("daniell", m = 50)) # A long moving average plot(kernel("daniell", c(5, 5))) # m=5 moving average of a m=5 moving average plot(kernel("daniell", c(5, 5, 5))) # a m=5 moving average of that! k = kernel("daniell", c(9, 9, 9)) smooth.spec <- spec.pgram(nino$SST, kernel = k, taper = 0) # Note how the confidence interval got much narrower spec.df <- data.frame(freq = smooth.spec$freq, `c(9,9,9)` = smooth.spec$spec) names(spec.df) <- c("freq", "c(9,9,9)") # Add other smooths k <- kernel("daniell", c(9, 9)) spec.df[, "c(9,9)"] <- spec.pgram(nino$SST, kernel = k, taper = 0, plot = FALSE)$spec k <- kernel("daniell", c(9)) spec.df[, "c(9)"] <- spec.pgram(nino$SST, kernel = k, taper = 0, plot = FALSE)$spec # melt from wide format into long format library(reshape2) spec.df <- melt(spec.df, id.vars = "freq", value.name = "spec", variable.name = "kernel") plot1 <- ggplot(data = subset(spec.df)) + geom_path(aes(x = freq, y = spec, color = kernel)) + scale_x_continuous("Period (years)", breaks = yrs.freqs, labels = yrs.labels) + scale_y_log10() plot2 <- ggplot(data = subset(spec.df)) + geom_path(aes(x = freq, y = spec, color = kernel)) + scale_x_log10("Period (years)", breaks = yrs.freqs, labels = yrs.labels) + scale_y_log10() grid.arrange(plot1, plot2) k = kernel("daniell", c(9, 9, 9)) smooth.spec <- spec.pgram(nino$SSA, kernel = k, taper = 0, plot = FALSE) # spec.df <- as.data.frame(smooth.spec) spec.df <- data.frame(freq = smooth.spec$freq, spec = smooth.spec$spec) ggplot(data = subset(spec.df)) + geom_line(aes(x = freq, y = spec)) + scale_x_continuous("Period (years)", breaks = yrs.freqs, labels = yrs.labels) + scale_y_continuous() k = kernel("daniell", c(9, 9)) smooth.spec <- spec.pgram(nino$SSA, kernel = k, taper = 0, plot = FALSE) spec.df <- data.frame(freq = smooth.spec$freq, `0%` = smooth.spec$spec) names(spec.df) <- c("freq", "0%") # Add other tapers spec.df[, "10%"] <- spec.pgram(nino$SSA, kernel = k, taper = 0.1, plot = FALSE)$spec spec.df[, "30%"] <- spec.pgram(nino$SSA, kernel = k, taper = 0.3, plot = FALSE)$spec spec.df <- melt(spec.df, id.vars = "freq", value.name = "spec", variable.name = "taper") plot1 <- ggplot(data = subset(spec.df)) + geom_path(aes(x = freq, y = spec, color = taper)) + scale_x_continuous("Period (years)", breaks = yrs.freqs, labels = yrs.labels) + scale_y_log10() plot2 <- ggplot(data = subset(spec.df)) + geom_path(aes(x = freq, y = spec, color = taper)) + scale_x_log10("Period (years)", breaks = yrs.freqs, labels = yrs.labels) + scale_y_log10() grid.arrange(plot1, plot2) k <- kernel("daniell", c(2)) spec.df[, "10%"] <- spec.pgram(nino$SSA, taper = 0.05)$spec library(multitaper) mt.spec <- spec.mtm(nino$SSA, nw = 16, k = 2 * 16 - 1, jackknife = TRUE, dtUnits = "month") # multitaper can resolve frequencies to about +/- NW/N Hz. i.e 16/1518 Hz # k is typically equal to 2NW - 1. Higher k is smoother mt.spec <- spec.mtm(nino$SST, nw = 16, k = 2 * 16 - 1, jackknife = TRUE, dT = 1/12, dtUnits = "year") library(dplR) wave.out <- morlet(y1 = nino$SST, x1 = nino$Year, p2 = 8, dj = 0.1, siglvl = 0.95) # p2=6 <=> estimate out to 2^8 = 256 months dj <=> controls the frequency # resolution hack the period estimate to be in years, not months wave.out$period <- wave.out$period/12 levs <- quantile(wave.out$Power, c(0, 0.25, 0.5, 0.75, 0.95, 1)) wavelet.plot(wave.out, wavelet.levels = levs, crn.ylim = c(22.5, 30)) wave.avg <- data.frame(power = apply(wave.out$Power, 2, mean), period = (wave.out$period)) plot(wave.avg$period, wave.avg$power, type = "l")