A cubic smoothing spline

Developed by Cook and Kairiukstis (1990), very popular and adjustable with the nyrs parameter.

spl32 <- ffcsaps(dat,nyrs=32)
spl60 <- ffcsaps(dat,nyrs=60)
spl128 <- ffcsaps(dat,nyrs=128)

par(tcl=0.5,mar=rep(3,4),mgp=c(1.1,0.1,0),xaxs="i")
plot(yrs,dat,type="l",xlab="Year",ylab="RWI",col="grey")
abline(h=1)
lines(yrs,spl32, col = "red", lwd = 2)
lines(yrs,spl60, col = "green", lwd = 2)
lines(yrs,spl128, col = "blue", lwd = 2)
axis(3);axis(4)

The Lowess Smoothing Spline

Piece wise linear polynomial fitting with smoothing width adjusted with f argument. Where f is the proportion of the data that is used in the smoothing.

f32 <- 32/nyrs
f32.lo <- lowess(x = yrs, y = dat, f = f32)
f60 <- 60/nyrs
f60.lo <- lowess(x = yrs, y = dat, f = f60)
f128 <- 128/nyrs
f128.lo <- lowess(x = yrs, y = dat, f = f128)

par(tcl=0.5,mar=rep(3,4),mgp=c(1.1,0.1,0),xaxs="i")
plot(yrs,dat,type="l",xlab="Year",ylab="RWI",col="grey")
abline(h=1)
lines(yrs,f32.lo$y, col = "red", lwd = 2)
lines(yrs,f60.lo$y, col = "green", lwd = 2)
lines(yrs,f128.lo$y, col = "blue", lwd = 2)
axis(3);axis(4)

Aggregating data

If you wanted to extract the median yearly values of a data set, say the December values of a time series you could use the filter function to extract or even aggregate something like the May through August values of a time series. The series below is the global ice extent. See Homework Section Below

Low Frequency variablity

Using RWI of tree growth data to remove low frequency variability in the data. Meaning the larger trend (low freq) where tree rings get smaller and smaller as a function of geometry as a tree ages. The first plot shows the decrease of ring with through age and the second plot shows the detrended data i.e. the RWI.

data(co021)
# the first tree-ring series
dat <- na.omit(co021[,"642244"])
n <- length(dat)
yrs <- 1:n

f128 <- 128/n
f128.lo <- lowess(x = yrs, y = dat, f = f128)
rwi <- dat/f128.lo$y

par(mfcol=c(2,1),mar=c(0.1,3,3,3),mgp=c(1.25,0.25,0),tcl=0.5)
plot(yrs,dat,type="l",ylab="mm",xlab="",col="grey",axes=F)
lines(yrs,f128.lo$y,lwd=2)
axis(2);axis(3);axis(4);box()
par(mar=c(3,3,0.1,3))
plot(yrs,rwi,type="l",ylab="RWI",xlab="Age",col="black",axes=F)
abline(h=1)
axis(2);axis(1);axis(4);box()

Fourier Transforms

To look at frequency data as opposed to time domain data. Break apart the data into as many sin and cos waves as necessary and tell R to look for the dominate or important frequencies in a time series that may be corrupted by noise. The following is a basic set of sine functions with a frequency of 5, 50, 100 and 250 Hz. We’ll then add some random noise.

n <- 1000
n.vec <- 1:n
wav5 <- 0.5 * sin(2 * pi / 5 * n.vec)
wav50 <- 0.75 * sin(2 * pi / 50 * n.vec)
wav100 <- sin(2 * pi / 100 * n.vec)
wav250 <- sin(2 * pi / 250 * n.vec)
plot(wav250,type="n")
lines(wav5,col="blue")
lines(wav50,col="green")
lines(wav100,col="red")
lines(wav250,col="purple")

noise <- rnorm(n)
wav <- wav250 + wav100 + wav50 + wav5 + noise 
plot(wav,type="l")

Using a Periodogram

To look at the break down of the major waves that appear when we break apart the set of sine waves we just made with the random noise. We use the spectrum function to do this. The lower the number or farther to left a peak is on the peridogram the lower its frequency is. remember \(frequency = \frac{1}{period}\) so for our largest period of 250, \(frequency = \frac{1}{250}=0.004\) and we see this as the most left peak in the following periodogram.

spectrum(wav,log="no")

Finally The Homework

First I will simply try to follow the first running average example just to make sure I have the format of my data in working order and well understood. I will use the northern ice extent data.

my.data<-readRDS("E:\\Flash Drive\\ESCI 597\\Data\\HW4\\ice.nx.rds")

yrs <- as.numeric(time(my.data))
dat <- my.data[1:458]
nyrs <- length(dat)

Ok so I finnnnaly got that to work and this doesn’t really show us anything particularly new that you can’t determine simply by visual inspection prior to the running a moving average.

ma32 <- filter(x=dat, filter=rep(x=1/32,times=32), sides=2)
ma60 <- filter(x=dat, filter=rep(x=1/60,times=60), sides=2)
ma128 <- filter(x=dat, filter=rep(x=1/128,times=128), sides=2)

par(tcl=0.5,mar=rep(3,4),mgp=c(1.1,0.1,0),xaxs="i")
plot(yrs,dat,type="l",xlab="Year",ylab="RWI",col="grey")
abline(h=1)
lines(yrs,ma32, col = "red", lwd = 2)
lines(yrs,ma60, col = "green", lwd = 2)
lines(yrs,ma128, col = "blue", lwd = 2)
axis(3);axis(4)

More Filtering

(From Above)

Because the smoothing and moving averages didn’t reveal much with the global sea ice extent, I will follow up and remove the winter three month maximums as a filtering and aggregating excercise. I downloaded data on the global ice extent from the world’s coolest Shiny app. The data below are monthly from Nov 1978 to December 2016. Let’s use the filter function to aggregate the data. We’ll get the average extent for winter months only (November, December and January).

Figuring out the alignment of the data for to make sure that I was actually getting the three month average of December (i.e. Nov, Dec, Jan) took for ever. I finnnnnallly got it figured out. Why do these things take soo long?

# extract the indices corresponding to each December, which will be the 
# average of Nov, Dec and Jan.
winter.index <-seq(from=2,to=length(ice.3.filt),by=12)
ice.winter <- ice.3.filt[winter.index]
# note the length of ice.summer
length(ice.winter)

## [1] 39

length(yrs)

## [1] 39

# compare
mean(ice.nx[1:3])

## [1] 13.92667

ice.winter[1]

## [1] 13.92667

par(mar=rep(4,4))

plot(yrs,ice.winter,type="b",xlab="Year",ylab=expression(1e6~km^2))
abline(lm(ice.winter~yrs),col="red")

Use the El Nino Data Instead

Lets try that with some more variable data like the el nino 3.4 data and see what it looks like.

The el nino data is much more nosiy and will work better for a lesson in spectral analysis. I will proceed with this data. Also, I made the plotted moving averages a tad smaller becasue they are very crowded with the origninal linewidth.

ma32 <- filter(x=dat, filter=rep(x=1/32,times=32), sides=2)
ma60 <- filter(x=dat, filter=rep(x=1/60,times=60), sides=2)
ma128 <- filter(x=dat, filter=rep(x=1/128,times=128), sides=2)

par(tcl=0.5,mar=rep(3,4),mgp=c(1.1,0.1,0),xaxs="i")
plot(yrs,dat,type="l",xlab="Year",ylab="RWI",col="grey")
abline(h=1)
lines(yrs,ma32, col = "red", lwd = 1)
lines(yrs,ma60, col = "green", lwd = 1)
lines(yrs,ma128, col = "blue", lwd = 1)
axis(3);axis(4)

I like the idea of piece wise functions so I will try a different spline instead of a filter creating a moving average and see if I can get a better look at what frequecies my be lurking in the el nino data. Using the lowess spline, I will try f=4, f=12 and f=32.

myts<-readRDS("E:\\Flash Drive\\ESCI 597\\Data\\HW4\\nino34.rds")

yrs <- as.numeric(time(myts))
dat <- myts[1:1762]
nyrs <- length(dat)

f4 <- 4/nyrs
f4.lo <- lowess(x = yrs, y = dat, f = f4)
f12 <- 12/nyrs
f12.lo <- lowess(x = yrs, y = dat, f = f12)
f32 <- 32/nyrs
f32.lo <- lowess(x = yrs, y = dat, f = f32)

par(tcl=0.5,mar=rep(3,4),mgp=c(1.1,0.1,0),xaxs="i")
plot(yrs,dat,type="l",xlab="Year",ylab="RWI",col="grey")
abline(h=1)
lines(yrs,f4.lo$y, col = "red", lwd = 2)
lines(yrs,f12.lo$y, col = "green", lwd = 2)
lines(yrs,f32.lo$y, col = "blue", lwd = 2)
axis(3);axis(4)

Still pretty variabel so increase the f values to f= 48, f= 60 and f=128.

Just for shits and giggles I will try one more spline because I want to see a smooth one. Polynomials are fun and cubic ones are even funner becuase they can provide information through the use of inflection points ina simple way so I will apply the ffcsaps method. Also I will bump the highest smoothing value up to say 256 just too see what that will look like. Note that 256 is in months so 256/12 would be just over 21 years.

Extracting the residuals for compairison.

The first set of plots is that of the el nino 3.4 index whith a 60 month loess smoothed spline (blue). The second set of plots is that of the el nino 3.4 index whith a 128 year loess smoothed spline (green). Both of the different loess smoothed time series are followed by a regression summary of thier residuals.

myts<-readRDS("E:\\Flash Drive\\ESCI 597\\Data\\HW4\\nino34.rds")

yrs <- as.numeric(time(myts))
dat <- myts[1:1762]
n <- length(dat)

f60 <- 60/n
f60.lo <- lowess(x = yrs, y = dat, f = f60)
nino60 <- dat/f60.lo$y

par(mfcol=c(2,1),mar=c(0.1,3,3,3),mgp=c(1.25,0.25,0),tcl=0.5)
plot(yrs,dat,type="l",ylab="Degrees Celsius",xlab="",col="grey",axes=F)
lines(yrs,f60.lo$y,lwd=2,col="blue")
axis(2);axis(3);axis(4);box()
par(mar=c(3,3,0.1,3))
plot(yrs,nino60,type="l",ylab="Detrended",xlab="Year",col="black",axes=F)
abline(h=1)
axis(2);axis(1);axis(4);box()

residuals<-lm(nino60~yrs)
summary(residuals)

## 
## Call:
## lm(formula = nino60 ~ yrs)
## 
## Residuals:
##       Min        1Q    Median        3Q       Max 
## -0.101629 -0.020472  0.001351  0.020880  0.076108 
## 
## Coefficients:
##               Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  1.029e+00  3.147e-02  32.695   <2e-16 ***
## yrs         -1.528e-05  1.619e-05  -0.944    0.345    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.02881 on 1760 degrees of freedom
## Multiple R-squared:  0.000506,   Adjusted R-squared:  -6.191e-05 
## F-statistic: 0.891 on 1 and 1760 DF,  p-value: 0.3453

The summary of the two different linear regression of the rsiduals show that the residuals of each different lowess smoothing spline are not drastically different as seen in the five number summary. With a 60 and 120 month loess spline we can somewhat see a 5 and 10 year ish year oscillation.

myts<-readRDS("E:\\Flash Drive\\ESCI 597\\Data\\HW4\\nino34.rds")

yrs <- as.numeric(time(myts))
dat <- myts[1:1762]
n <- length(dat)

f120 <- 120/n
f120.lo <- lowess(x = yrs, y = dat, f = f120)
nino120 <- dat/f120.lo$y

par(mfcol=c(2,1),mar=c(0.1,3,3,3),mgp=c(1.25,0.25,0),tcl=0.5)
plot(yrs,dat,type="l",ylab="Degrees Celsius",xlab="",col="grey",axes=F)
lines(yrs,f120.lo$y,lwd=2,col="green")
axis(2);axis(3);axis(4);box()
par(mar=c(3,3,0.1,3))
plot(yrs,nino120,type="l",ylab="Detrended",xlab="Year",col="black",axes=F)
abline(h=1)
axis(2);axis(1);axis(4);box()

residuals<-lm(nino120~yrs)
summary(residuals)

## 
## Call:
## lm(formula = nino120 ~ yrs)
## 
## Residuals:
##       Min        1Q    Median        3Q       Max 
## -0.097765 -0.022695  0.000709  0.021758  0.086321 
## 
## Coefficients:
##               Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  1.029e+00  3.417e-02  30.129   <2e-16 ***
## yrs         -1.526e-05  1.758e-05  -0.868    0.386    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.03127 on 1760 degrees of freedom
## Multiple R-squared:  0.0004279,  Adjusted R-squared:  -0.0001401 
## F-statistic: 0.7534 on 1 and 1760 DF,  p-value: 0.3855

Third I will filter a methane ch4.grim-1.rds data such that it will be transformed from monthly to yearly data. It is not a very long data set but will serve the purpose of creating a annual data set from a monthly ts.

## [1] 464

Filtering to create an annual data set

From a monthly data set. This particular data set starts in April. Therefor six months later would be September and I centered a 12 month average on September and used sequence of numbers starding ot the 6th position in this data set to extract the mean ch4 value at every September and then plot that data.

yrs <- 1979:2016
# This is 3-month filter 
ch4.12.filt <- filter(x = ch4,filter = rep(x=1/12,times=12),sides=2)
# Note the NA
head(cbind(ch4,ch4.12.filt))

##           ch4 ch4.12.filt
## [1,] 1479.000          NA
## [2,] 1485.033          NA
## [3,] 1491.067          NA
## [4,] 1497.100          NA
## [5,] 1496.814          NA
## [6,] 1496.529    1493.447

# extract the indices corresponding to each December, which will be the average of Dec - Dec.
ch4.index <-seq(from=6,to=length(ice.3.filt),by=12)
annual.ch4 <- ch4.12.filt[ch4.index]
# note the length of ice.summer
length(annual.ch4)

## [1] 38

length(yrs)

## [1] 38

# compare
mean(ch4[1:12])

## [1] 1493.447

annual.ch4[1]

## [1] 1493.447

par(mar=rep(4,4))

plot(yrs,annual.ch4,type="b",xlab="Year",ylab="PPB")
abline(lm(annual.ch4~yrs),col="red")

Spectral analysis

Orbital forcing and Milankovitch cycles

To Be Finished very early in the morning:(

jul65N <- readRDS("E:\\Flash Drive\\ESCI 597\\Data\\HW4\\jul65N.rds")
plot(jul65N$kya,jul65N$W.per.m2,type="l",
     xlab="Thousands of years ago",
     ylab=expression(W/m^2),main="July insolation at 65N")

ss.spec <- spectrum(jul65N, span=2, plot= FALSE)
plot(ss.spec,log="no", type="h", xlab="Frequency (cycles / yr)", ylab="Spectral density",main="Annual solar insolation")

plot(ss.spec, log="no", type="h", xlab="Frequency (cycles / yr)", ylab="Spectral density",main="AAnnual solar insolation",xlim=c(0,0.08))
abline(v=0.0245)
abline(v=0.042)
abline(v=0.045)
abline(v=0.053)

Sooooo if I understand this correctly the dominate frequencies that are occuring in the annual insolation data are as follows:

*1/0.0245=40.82

*1/0.04=25

*1/0.045=22.2

*1/0.053=18.87

These values correspond nicely to the Milankovtich cycles that are related to the earths tilt (41K year cycles), axial precession (21K year cycles) and apisidal precession (26K year cycles) nicely! That is THE NUTS!!! and super cool. If I had time I would love to look more into the mechanics of all the orbital physics that underpin this.

HW4_Klinger

Moving Average filter

Smoothing Splines