--- title: " " output: html_document: default pdf_document: default ---

# Theory and Applications

### SSP Shen and RCJ Somerville

Version 1.0 released in July 2019 and coded by Dr. Samuel Shen, Distinguished Professor
San Diego State University, California, USA
Email:

Version 2.0 compiled and updated by Momtaza Sayd
San Diego State University May 2021
Version 3.0 compiled and updated by Joaquin Stawsky
San Diego State University June 2022

## Chapter 4: Climate Data Matrices and Linear Algebra ### Matrix algebra ```{r} #Matrix subtraction matrix(c(1,1,1,-1), nrow = 2) - matrix(c(1,3,2,0), nrow = 2) #Matrix multiplication matrix(c(1,1,1,-1), nrow = 2) %*% matrix(c(1,2,3,4), nrow = 2) matrix(c(1,2,3,4), nrow = 2) %*% matrix(c(1,1,1,-1), nrow = 2) #Find the inverse of a matrix solve(matrix(c(1,1,1,-1), nrow = 2)) #Verify the inverse matrix matrix(c(0.5,0.5,0.5,-0.5), nrow = 2) %*% matrix(c(1,1,1,-1), nrow = 2) #Solve a group of linear equations solve(matrix(c(1,1,1,-1),nrow = 2),c(20,4)) #This is the result x1 = 12, and x2 = 8. #Compute a covariance matrix A <- matrix(c(1,-1,2,0,3,1),nrow = 2) A covm <- (1/(dim(A)))*A%*%t(A) covm #is the covariance matrix. #Matrix multiplication of a vector changes the vector direction u <- c(1,-1) v <- covm%*%u v #u and v are in different directions. #Matrix mutiplication of an eigenvector does change direction ew <- eigen(covm) eigen(covm)\$values eigen(covm)\$vectors #Verify the eigenvectors and eigenvalues covm%*%ew\$vectors[,1]/ew\$values #This is the first eigenvector #SVD: Decomposition and recovery of a space-time data matrix #Develop a 2-by-3 space-time data matrix for SVD A <- matrix(c(1,-1,2,0,3,1),nrow = 2) A #Perform SVD calculation msvd <- svd(A) msvd msvd\$d msvd\$u msvd\$v #One can verify that A = UDV', where V' is transpose of V. verim <- msvd\$u%*%diag(msvd\$d)%*%t(msvd\$v) verim round(verim) #This is the original data matrix A #Covariance matrix and its eigen-problem covm <- (1/(dim(A)))*A%*%t(A) eigcov <- eigen(covm) eigcov\$values eigcov\$vectors #Eigenvalues from SVD and from a covariance matrix ((msvd\$d)^2)/(dim(A)) eigcov\$values ```

### Plot Fig. 4.2 ```{r} #SOI and the Standardized SLP data at Darwin and Tahiti setwd("~/sshen/climmath") #setEPS() #postscript("fig0402a.eps", height = 4.9, width = 7) #par(mar = c(4.0,4.2,1.5,0.5)) Pta <- read.table("data/PSTANDtahiti.txt", header = F) # Remove the first column that is the year ptamon <- Pta[,seq(2,13)] #Convert the matrix into a vector according to mon: Jan 1951, Feb 1951, ..., Dec 2015 ptamonv <- c(t(ptamon)) #Generate time ticks from Jan 1951 to Dec 2015 xtime <- seq(1951, 2016-1/12, 1/12) # Plot the Tahiti standardized SLP anomalies plot(xtime, ptamonv,type = "l",xlab = "Year", ylab = "Pressure Anomalies", cex.lab = 1.2, cex.axis = 1.2, main = "Standardized Tahiti SLP Anomalies", col = "red", xlim = range(xtime), ylim = c(-4,4)) text(1952,3.5, "(a)", cex = 1.2) # dev.off() ```

### Do the same for Darwin SLP ```{r} setwd("~/sshen/climmath/") #setEPS() #postscript("fig0402b.eps", height = 4.9, width = 7) #par(mar = c(4.0,4.2,1.5,0.5)) Pda <- read.table("data/PSTANDdarwin.txt", header = F) pdamon <- Pda[,seq(2,13)] pdamonv <- c(t(pdamon)) plot(xtime, pdamonv,type = "l",cex.lab = 1.2, cex.axis = 1.2, xlab = "Year",ylab = "Pressure", main = "Standardized Darwin SLP Anomalies", col = "blue", xlim = range(xtime), ylim = c(-4,4)) text(1952,3.5, "(b)", cex = 1.2) # dev.off() ```

### Plot the SOI index ```{r} #setEPS() #postscript("fig0402c.eps", height = 4.9, width = 7) #par(mar = c(4.0,4.2,2.0,0.5)) SOI <- ptamonv-pdamonv plot(xtime, SOI ,type = "l", cex.lab = 1.2, cex.axis = 1.2, xlab = "Year",ylab = "SOI index", col = "black",xlim = range(xtime), ylim = c(-6,6), lwd = 1) #Add ticks on top edge of the plot box axis (3, at = seq(1951,2015,4), labels = seq(1951,2015,4)) #Add ticks on the right edge of the plot box axis (4, at = seq(-4,4,2), labels = seq(-4,4,2)) lines(xtime, rep(0,length(xtime))) text(1985, 5, "SOI = Tahiti Pressure - Darwin Pressure", cex = 1.2) text(1952,5.6, "(c)", cex = 1.2) #abline(lm(SOI ~ xtime), col = "red", lwd = 2) # dev.off() ```

### CSOI and AMO time series comparison ```{r} setwd("~/sshen/climmath/") # setEPS() #postscript("fig0402d.eps", height = 4.9, width = 7) par(mar = c(4,4.2,1.6,4)) cnegsoi <- -cumsum(ptamonv-pdamonv) #par(mgp = c(2,2,4,0)) plot(xtime, cnegsoi,type = "l", cex.lab = 1.2, cex.axis = 1.2, xlab = "Year",ylab = "Negative CSOI index", main = "CSOI and AMO Index Comparison", col = "purple",xlim = range(xtime), ylim = range(cnegsoi), lwd = 1.5) legend(1960,15, col = c("purple"),lty = 1,lwd = 2.0, legend = c("CSOI"),bty = "n",text.font = 2,cex = 1.2) text(1951,12, "(d)", cex = 1.2) #AMO data and plot amodat <- read.table("data/AMO1951-2015.txt", header = FALSE) amots <- as.vector(t(amodat[,2:13])) par(new = TRUE) plot(xtime, amots,type = "l",col = "darkgreen", cex.lab = 1.2, cex.axis = 1.2, lwd = 1.5,axes = FALSE,xlab = "",ylab = "") legend(1960,0.45, col = c("darkgreen"),lty = 1,lwd = 2.0, legend = c("AMO index"),bty = "n",text.font = 2,cex = 1.2) #Suppress the axes and assign the y-axis to side 4 axis(4, cex.axis = 1.3) mtext("AMO index",side = 4,line = 3, cex = 1.2) # dev.off() ```

### Compute the weighted SOI by SVD ```{r} #Time-space data format ptada <- cbind(ptamonv,pdamonv) #Space-time data format ptada <- t(ptada) #Implement SVD for the sapce-time data svdptd <- svd(ptada) recontd <- svdptd\$u%*%diag(svdptd\$d[1:2])%*%t(svdptd\$v) recontd <- ptada #Two EOFs from the space-time data U <- svdptd\$u U #These are the weights for the WSOI svdptd\$d D <- diag(svdptd\$d) D #PCs: The principal components V <- svdptd\$v V ```

### Plot Fig. 4.3 ```{r} U <- svdptd\$u #The matrix of the EOF spatial patterns V <- svdptd\$v #The matrix of the PC temporal patterns D <- diag(svdptd\$d) #The diagonal matrix of standard deviations recontada <- U%*%D%*%t(V) #Verify the SVD decomposition svd_error <- ptada - recontada max(svd_error) min(svd_error) ```

### Plot WSOI1 ```{r} #setEPS() #postscript("fig0403a.eps", height = 4.9, width = 7) par(mar = c(4,4.2,2.0,0.5)) xtime <- seq(1951, 2016-1/12, 1/12) wsoi1 <- D[1,1]*t(V)[1,] plot(xtime, wsoi1,type = "l", cex.lab = 1.2, cex.axis = 1.2, xlab = "Year",ylab = "WSOI1", col = "black",xlim = range(xtime), ylim = c(-4,4)) axis (3, at = seq(1951,2015,4), labels = seq(1951,2015,4)) text(1985,-3.7,"Weighted SOI1 Index", cex = 1.2) text(1952,3.5, "(a)", cex = 1.2) # dev.off() ```

### Plot WSOI2 ```{r} #setEPS() #postscript("fig0403b.eps", height = 4.9, width = 7) par(mar = c(4,4.2,2.0,0.5)) xtime <- seq(1951, 2016-1/12, 1/12) wsoi2 <- D[2,2]*t(V)[2,] plot(xtime, wsoi2,type = "l", cex.lab = 1.3, cex.axis = 1.3, xlab = "Year",ylab = "WSOI2", col = "blue",xlim = range(xtime), ylim = c(-4,4)) axis (3, at = seq(1951,2015,4), labels = seq(1951,2015,4)) text(1985,-3.7,"Weighted SOI2 Index", cex = 1.5) text(1952,3.5, "(b)", cex = 1.3) # dev.off() ```

### Cumulative WSOIs ```{r} setwd("~/sshen/climmath/") #setEPS() #postscript("fig0403c.eps", height = 4.9, width = 7) par(mar = c(4,4.2,2.0,4.2)) cwsoi1 <- cumsum(wsoi1) plot(xtime, cwsoi1,type = "l",cex.lab = 1.2, cex.axis = 1.2, xlab = "Year",ylab = "Weighted SOI1", col = "black",lwd = 3, ylim = c(-200,50), main = "Comparison between CWSOI1 and the Smoothed AMO Index") text(1951,40, "(c)", cex = 1.3) #axis (3, at = seq(1951,2015,4), labels = seq(1951,2015,4)) legend(1970,20, col = c("black"),lty = 1,lwd = 3.0, legend = c("CWSOI1"),bty = "n",text.font = 2,cex = 1.2) #Superimpose CSOI time series on this CWSOI1 cnegsoi <- -cumsum(ptamonv-pdamonv) lines(xtime, cnegsoi,type = "l",col = "purple", lwd = 3.0) legend(1970,2, col = c("purple"),lty = 1,lwd = 3.0, legend = c("CSOI"),bty = "n",text.font = 2,cex = 1.2) #24-month ahead moving average of the monthly AMO index amodat <- read.table("data/AMO1951-2015.txt", header = FALSE) amots <- as.vector(t(amodat[,2:13])) #install.packages("TTR") library("TTR") amomv <- SMA(amots,n = 24, fill = NA) #Average of the previous n points par(new = TRUE) xtime <- seq(1951,2015,len = 780) plot(xtime, amomv[37:816],type = "l",col = "darkgreen", lwd = 3,axes = FALSE,xlab = "",ylab = "", cex.lab = 1.2, cex.axis = 1.2, ylim = c(-0.5, 0.5)) legend(1970,0.23, col = c("darkgreen"),lty = 1,lwd = 3.0, legend = c("AMO index"),bty = "n",text.font = 2,cex = 1.2) #Suppress the axes and assign the y-axis to side 4 axis(4, cex.axis = 1.2, col.axis = "darkgreen") mtext("Smoothed AMO index", side = 4,line = 3, cex = 1.2, col = "darkgreen") # dev.off() ```

### Plot cumulative WSOI2: CWSOI2 ```{r} #setEPS() #postscript("fig0403d.eps", height = 4.9, width = 7) par(mar = c(4,4.2,2.0,0.5)) cwsoi1 <- cumsum(wsoi1) wsoi2 <- D[2,2]*t(V)[2,] cwsoi2 <- cumsum(wsoi2) plot.new() # par(mar = c(4,4,3,4)) plot(xtime, cwsoi2,type = "l",xlab = "Year",ylab = "CWSOI2", cex.lab = 1.2, cex.axis = 1.2, col = "blue",lwd = 3, ylim = c(-200,50), main = "CWSOI2 Index: The Cumulative PC2") text(1951,40, "(d)", cex = 1.3) # dev.off() ```

### Plot Fig. 4.4 ```{r} #Display the two ENSO modes on a world map library(maps) library(mapdata) plot.new() par(mfrow = c(2,1)) par(mar = c(0,0,0,0)) #Zero space between (a) and (b) map(database = "world2Hires",ylim = c(-70,70), mar = c(0,0,0,0)) grid(nx = 12,ny = 6) points(231, -18,pch = 16,cex = 2, col = "red") text(231, -30, "Tahiti 0.61", col = "red") points(131, -12,pch = 16,cex = 2.6, col = "blue") text(131, -24, "Darwin -0.79", col = "blue") axis(2, at = seq(-70,70,20), col.axis = "black", tck = -0.05, las = 2, line = -0.9,lwd = 0) axis(1, at = seq(0,360,60), col.axis = "black",tck = -0.05, las = 1, line = -0.9,lwd = 0) text(180,30, "El Nino Southern Oscillation Mode 1",col = "purple",cex = 1.3) text(10,-60,"(a)", cex = 1.4) box() # par(mar = c(0,0,0,0)) #Plot mode 2 map(database = "world2Hires", ylim = c(-70,70), mar = c(0,0,0,0)) grid(nx = 12,ny = 6) points(231, -18,pch = 16,cex = 2.6, col = "red") text(231, -30, "Tahiti 0.79", col = "red") points(131, -12,pch = 16,cex = 2, col = "red") text(131, -24, "Darwin 0.61", col = "red") text(180,30, "El Nino Southern Oscillation Mode 2",col = "purple",cex = 1.3) axis(2, at = seq(-70,70,20), col.axis = "black", tck = -0.05, las = 2, line = -0.9,lwd = 0) axis(1, at = seq(0,360,60), col.axis = "black",tck = -0.05, las = 1, line = -0.9,lwd = 0) text(10,-60,"(b)", cex = 1.4) box() ```

### An example of multivariate regression ```{r} x1 <- c(1,2,3) #Given the coordinates of the 3 points x2 <- c(2,1,3) y <- c(-1,2,1) df <- data.frame(x1,x2,y) #Put data into the data.frame format fit <- lm(y ~ x1 + x2, data = df) fit #Show the regression results 1.667*x1-1.333*x2 #Verify that 3 points determining a plane ```

### Multivariate regression with four data points ```{r} u <- c(1,2,3,1) v <- c(2,4,3,-1) w <- c(1,-2,3,4) mydata <- data.frame(u,v,w) myfit <- lm(w ~ u + v, data = mydata) summary(myfit) #Show the result ```

### Another example of multivariate regregression ```{r} dat <- matrix(rnorm(40),nrow = 10, dimnames = list(c(letters[1:10]), c(LETTERS[23:26]))) fdat <- data.frame(dat) fit <- lm(Z ~ W + X + Y, data = fdat) summary(fit) ```