##########################################
# FILE: EVModelingMAMps.txt - SPRING
# AUTHOR: Elizabeth Shamseldin
# DATE: 3/7/06
# PURPOSE: Create point-level yr-return estimates modeled from grid-level yr-returns
# using parameter estimates from
# http://www.stat.unc.edu/postscript/rs/ext1/ for NCAR point vs grid-cell models
# Program CLEANS data for |xi|>1, logpsi, mu, grid not matching point, non-convergence
# 
# NOTE: On Elizabeth's Mac, X11 has to be RUNNING in order to produce png plots
#
# MODIFIED: 3/14/06
# 3/9: Automated program to export models and plots to PNG files (not pdf)
# 3/14: Changed pdf to png files
# 3/14: set working directories PNG and ModelParams
# 3/14: Included Residual degrees of freedom in saved Model output
# 3/14: Cleaned up output/file-naming convention and folder system
#	 (swapped p.returnVSg.return for 50Ret and 100Ret folders)
#
# 3/17: Modified desktop program to run on Fisher
#		Changed all setwd() references
#		First run on "old" data, then will run on new param estimates
#
# 3/20: "Old" param estimates moved to EVdataOLD
#		"NEW" param estimates downloaded from http://www.stat.unc.edu/postscript/rs/ext1/file3/
#		"NEW" estimates in EVdata
##########################################



library(fields)
library(evd)
library(ismev)

setwd("/Users/shamseld/Desktop/NCAR/")



#Set Confidence level z-value for plots
CI.z <- 1.96
#Set Header for Model Output
param.head <- cbind("1 intercept", "X1",  "aic",  "deviance", "numObs")
param.headE <- cbind("1 intercept", "X1",  "X2(Elev)", "aic",  "deviance", "numObs")
param.head50100 <- cbind("1 intercept50", "X50-50", "X50-100", "intercept100", "X100-50", "X100-100", "numObs")
param.head50100E1 <- cbind("1 intercept50", "X50-50", "X50-100", "X50(Elev)", 
	"intercept100", "X100-50", "X100-100", "X100(Elev)", "numObs")
param.head50100E2 <- cbind("1 intercept50", "X50-50", "X50-100", "X50(Elev50)", "X50(Elev100)", 
	"intercept100", "X100-50", "X100-100", "X100(Elev50)", "X50(Elev100)",  "numObs")


#Header
head <- cbind("num_station", "num_obs", "thresh", "lat", "lon", "elev", "conv", "mu_mle", "logpsi_mle", "xi_mle", "mu_se", "logpsi_se", "xi_se", "mulogpsi_corr", "muxi_corr", "logpsixi_corr")

#Read in grid data.  sep=""
g.DATA <- read.table(file="EVdata/GDJF19501999109511.txt", sep="", header=FALSE,  col.names=head) 
p.DATA <- read.table(file="EVdata/PDJF19501999109511.txt", sep="", header=FALSE,  col.names=head) 

StartTime <- Sys.time()
# Length: 4342
g.gridnum <- g.DATA[,1]
g.lat <- g.DATA[,4]
g.latmin <- g.lat-1.25
g.latmax <- g.lat+1.25
g.lon <- g.DATA[,5]
g.lonmin <- g.lon-1.25
g.lonmax <- g.lon+1.25
g.mu <- g.DATA[,8]
g.logpsi <- g.DATA[,9]
g.xi <- g.DATA[,10]
g.conv <- g.DATA[,7]

g.mu_se <- g.DATA[,11]
g.logpsi_se <- g.DATA[,12]
g.xi_se <- g.DATA[,13]
g.mulogpsi_corr <- g.DATA[,14]
g.muxi_corr <- g.DATA[,15]
g.logpsixi_corr <- g.DATA[,16]


p.lat <- p.DATA[,4]
p.lon <- p.DATA[,5]
p.elev <- p.DATA[,6]
p.mu <- p.DATA[,8]
p.logpsi <- p.DATA[,9]
p.xi <- p.DATA[,10]
p.gridnum <- matrix(0,length(p.DATA[,1]),1)
p.grid.lat <- matrix(0,length(p.DATA[,1]),1)
p.grid.lon <- matrix(0,length(p.DATA[,1]),1)

p.grid.mu <- matrix(0,length(p.DATA[,1]),1)
p.grid.logpsi <- matrix(0,length(p.DATA[,1]),1)
p.grid.xi <- matrix(0,length(p.DATA[,1]),1)
p.conv <- p.DATA[,7]
p.grid.conv <- matrix(0, length(p.DATA[,1]), 1)

p.grid.mu_se <- matrix(0,length(p.DATA[,1]),1)
p.grid.logpsi_se <- matrix(0,length(p.DATA[,1]),1)
p.grid.xi_se <- matrix(0,length(p.DATA[,1]),1)
p.grid.mulogpsi_corr <- matrix(0,length(p.DATA[,1]),1)
p.grid.muxi_corr <- matrix(0,length(p.DATA[,1]),1)
p.grid.logpsixi_corr <- matrix(0,length(p.DATA[,1]),1)

#Assign grids to point locations via lat and lon
for(i in 1:length(g.DATA[,1]))
	{
	for(j in 1:length(p.DATA[,1]))
		{
			if(p.lat[j] >= g.latmin[i] && p.lat[j] <= g.latmax[i] && p.lon[j] >= g.lonmin[i] && p.lon[j] <= g.lonmax[i])
			{p.gridnum[j] <- g.gridnum[i]}

		}#end j.p for
	}#end i.g for


#Assign mle values to point matrix
for(i in 1:length(p.DATA[,1]))
	{
		p.grid.lat[i] <- g.lat[(p.gridnum[i])]
		p.grid.lon[i] <- g.lon[(p.gridnum[i])]
		p.grid.mu[i] <- g.mu[(p.gridnum[i])]
		p.grid.logpsi[i] <- g.logpsi[(p.gridnum[i])]
		p.grid.xi[i] <- g.xi[(p.gridnum[i])]
		p.grid.conv[i] <- g.conv[(p.gridnum[i])]
		p.grid.mu_se[i] <- g.mu_se[(p.gridnum[i])]
		p.grid.logpsi_se[i] <- g.logpsi_se[(p.gridnum[i])]
		p.grid.xi_se[i] <- g.xi_se[(p.gridnum[i])]
		p.grid.mulogpsi_corr[i] <- g.mulogpsi_corr[(p.gridnum[i])]
		p.grid.muxi_corr[i] <- g.muxi_corr[(p.gridnum[i])]
		p.grid.logpsixi_corr[i] <- g.logpsixi_corr[(p.gridnum[i])]

	}


########################################################################################
########################################################################################
########################################################################################
#CLEANING DATA
# Program CLEANS data for |xi|>1, logpsi, mu, grid not matching point, non-convergence

Sys.time()  #Takes about 15 minutes

l <- length(p.DATA[,1])

p.TEMP <- cbind(p.DATA[1:l,1:16], p.gridnum[1:l], p.grid.conv[1:l], p.grid.lat[1:l], p.grid.lon[1:l], p.grid.mu[1:l], 
		p.grid.logpsi[1:l], p.grid.xi[1:l], p.grid.mu_se[1:l], p.grid.logpsi_se[1:l], p.grid.xi_se[1:l], 
		p.grid.mulogpsi_corr[1:l], p.grid.muxi_corr[1:l], p.grid.logpsixi_corr[1:l])


p.TEMP2 <- NULL

#for(i in 1:length(p.DATA[,1])){
for(i in 1:l){

	if( p.grid.xi[i] < 1 && p.grid.xi[i] > -1 && p.xi[i] < 1 && p.xi[i]  > -1 &&
		p.grid.logpsi[i] > 2 && p.grid.logpsi[i] < 6 && p.logpsi[i] > 2 && p.logpsi[i] < 6 &&
		p.grid.conv[i] == 0 && p.conv[i] == 0 ){

		p.TEMP2 <- rbind(p.TEMP2, p.TEMP[i,])
	}#end if

}#end for

Sys.time()

#setwd("/Users/shamseld/Desktop/NCAR/PointVSGrid/WINTER")
write.table(p.TEMP2, file="PointVSGrid/WINTER/PDJFCl.csv")

########################################################################################
########################################################################################
########################################################################################
# Length: 4339
pc.lat <- p.TEMP2[,4]
pc.lon <- p.TEMP2[,5]
pc.elev <- p.TEMP2[,6]
pc.mu <- p.TEMP2[,8]
pc.logpsi <- p.TEMP2[,9]
pc.xi <- p.TEMP2[,10]
pc.gridnum <- p.TEMP2[,17]
pc.grid.lat <- p.TEMP2[,19]
pc.grid.lon <- p.TEMP2[,20]

pc.grid.mu <- p.TEMP2[,21]
pc.grid.logpsi <- p.TEMP2[,22]
pc.grid.xi <- p.TEMP2[,23]
pc.conv <- p.TEMP2[,7]
pc.grid.conv <- p.TEMP2[,18]

pc.grid.mu_se <- p.TEMP2[,24]
pc.grid.logpsi_se <- p.TEMP2[,25]
pc.grid.xi_se <- p.TEMP2[,26]
pc.grid.mulogpsi_corr <- p.TEMP2[,27]
pc.grid.muxi_corr <- p.TEMP2[,28]
pc.grid.logpsixi_corr <- p.TEMP2[,29]


pc.return100 <- matrix(0, length(p.TEMP2[,1]), 1)
pc.grid.return100 <- matrix(0, length(p.TEMP2[,1]), 1)
n.return <- 100

#Set return level based on xi=0, xi!=0
for(i in 1:length(p.TEMP2[,1]))
	{
		if(pc.xi[i] != 0) 
			{ pc.return100[i] <- (pc.mu[i] + exp(pc.logpsi[i])*((n.return^pc.xi[i] - 1)/pc.xi[i])) }
		if(pc.grid.xi[i] != 0)
 			{ pc.grid.return100[i] <- (pc.grid.mu[i] + exp(pc.grid.logpsi[i])*((n.return^pc.grid.xi[i] - 1)/pc.grid.xi[i])) }
		
		if(pc.xi[i] == 0)
			{ pc.return100[i] <- (pc.mu[i] + exp(pc.logpsi[i])*log(n.return)) }
 		if(pc.grid.xi[i] == 0)
 			{ pc.grid.return100[i] <- (pc.grid.mu[i] + exp(pc.grid.logpsi[i])*log(n.return)) }
 		
	}


########################################################################################
########################################################################################
#setwd(wdPNG)
png(filename="PointVSGrid/WINTER/Models/GLM/100Ret/PNG/100RetDJFPtVSGrid.png")
plot(pc.grid.return100, pc.return100, xlab="Grid 100-yr Return", ylab="Point 100-yr Return", cex=0.5, col="green", 
main="Point (Green) vs Grid (Black) 100-Yr Return DJF")
for(i in 1:length(p.TEMP2[,1]))
	{
		points(pc.grid.return100[i], pc.grid.return100[i], col = "black", cex=0.1)
	} 
dev.off()
########################################################################################


return100c.glm <- glm(pc.return100 ~ pc.grid.return100)   
pred100 <- predict(return100c.glm, se.fit=TRUE)
out100 <- c(return100c.glm$coefficients, return100c.glm$aic, return100c.glm$deviance, length(return100c.glm$residuals))
#setwd(wdMP)
write.table(t(out100), file="PointVSGrid/WINTER/Models/GLM/100Ret/ModelParams/100RetDJFModel.csv", col.names=param.head)

#setwd(wdPNG)
png(filename="PointVSGrid/WINTER/Models/GLM/100Ret/PNG/100RetResidDJF.png", pointsize=6)
plot(pc.grid.return100, return100c.glm$residuals, main="100-Yr Ret Residuals - DJF", cex=0.5)
dev.off()


png(filename="PointVSGrid/WINTER/Models/GLM/100Ret/PNG/100RetStuResidDJF.png")
plot(pc.grid.return100, rstudent(return100c.glm), main="100-Yr Ret Student Residuals - DJF", cex=0.5)
dev.off()

png(filename="PointVSGrid/WINTER/Models/GLM/100Ret/PNG/100RetFitDJF.png")
plot(pc.grid.return100, pc.return100, main="100-Yr Ret Fitted Values DJF", cex=0.5, col="green")
for(i in 1:length(p.TEMP2[,1]))
	{
		points(pc.grid.return100[i], return100c.glm$fitted.values[i], cex=0.4)
		points(pc.grid.return100[i], (return100c.glm$fitted.values[i]+CI.z*pred100$se.fit[i]), cex = 0.4, col="red")
		points(pc.grid.return100[i], (return100c.glm$fitted.values[i]-CI.z*pred100$se.fit[i]), cex = 0.4, col="red")
	}
dev.off() 



###########################################################################################
###########################################################################################
###########################################################################################

return100LOGc.glm <- glm(log(pc.return100) ~ pc.grid.return100)  
pred100LOG <- predict(return100LOGc.glm, se.fit=TRUE)
out100LOG <- c(return100LOGc.glm$coefficients, return100LOGc.glm$aic, return100LOGc.glm$deviance, length(return100LOGc.glm$residuals))
#setwd(wdMP)
write.table(t(out100LOG), file="PointVSGrid/WINTER/Models/GLM/100Ret/ModelParams/100RetDJFLOGPtModel.csv", col.names=param.head)

#setwd(wdPNG)
png(filename="PointVSGrid/WINTER/Models/GLM/100Ret/PNG/100RetLogPtResidDJF.png")
plot(pc.grid.return100, return100LOGc.glm$residuals, main="100-Yr Ret LogPt Resid - DJF", cex=0.4)
dev.off() 

png(filename="PointVSGrid/WINTER/Models/GLM/100Ret/PNG/100RetLogPtStuResidDJF.png")
plot(pc.grid.return100, rstudent(return100LOGc.glm), main="100-Yr Ret LogPt StuResid - DJF", cex=0.4)
dev.off() 

png(file="PointVSGrid/WINTER/Models/GLM/100Ret/PNG/100RetLogPtFitDJF.png")
plot(pc.grid.return100, log(pc.return100), main="100-Yr Return  LogPt Fitted Values - DJF", cex=0.4, col="green")
for(i in 1:length(p.TEMP2[,1]))
	{
		points(pc.grid.return100[i], return100LOGc.glm$fitted.values[i], cex=0.5)
		points(pc.grid.return100[i], (return100LOGc.glm$fitted.values[i]+CI.z*pred100LOG$se.fit[i]), cex = 0.4, col="red")
		points(pc.grid.return100[i], (return100LOGc.glm$fitted.values[i]-CI.z*pred100LOG$se.fit[i]), cex = 0.4, col="red")
	}
dev.off() 
########################################################################
########################################################################
########################################################################
#Clean studentized residuals magnitude greater than 2
#Weight based on studentized residuals
pc.return.weight <- matrix(1, length(p.TEMP2[,1]), 1)

for(i in 1:length(p.TEMP2[,1]))  #This takes ~10 minutes to run
	{
		if(rstudent(return100LOGc.glm)[i] > 2 || rstudent(return100LOGc.glm)[i] < -2)
		{pc.return.weight[i] <- 0}

	}

pc.grid.return100W <- NULL
pc.return100W <- NULL
for(i in 1:length(p.TEMP2[,1]))  
	{
		if(pc.return.weight[i] != 0)
		{pc.grid.return100W <- rbind(pc.grid.return100W, pc.grid.return100[i]);
		pc.return100W <- rbind(pc.return100W, pc.return100[i])}

	}

#setwd("/Users/shamseld/Desktop/NCAR/PointVSGrid/WINTER")
write.table(cbind(pc.grid.return100W, pc.return100W), file="PointVSGrid/WINTER/PDJF100W.csv")

return100LOGcW.glm <- glm(log(pc.return100W) ~ pc.grid.return100W) 
pred100LOGW <- predict(return100LOGcW.glm, se.fit=TRUE, )
out100LOGW <- c(return100LOGcW.glm$coefficients, return100LOGcW.glm$aic, return100LOGcW.glm$deviance,
length(return100LOGcW.glm$residuals))
#setwd(wdMP)
write.table(t(out100LOGW), file="PointVSGrid/WINTER/Models/GLM/100Ret/ModelParams/100RetDJFLOGPtModelW.csv", col.names=param.head)

#setwd(wdPNG)
png(file="PointVSGrid/WINTER/Models/GLM/100Ret/PNG/100RetLogPtStuResidDJFW.png")
plot(pc.grid.return100W, rstudent(return100LOGcW.glm), main="100-Yr Ret LogPt (W) Student Resid - DJF", cex=0.4)
dev.off() 

png(file="PointVSGrid/WINTER/Models/GLM/100Ret/PNG/100RetLogPtFitDJFW.png")
plot(pc.grid.return100W, log(pc.return100W), main="100-Yr Ret LogPt (W) Fitted Values - DJF", cex=0.4, col="green")
for(i in 1:length(p.TEMP2[,1]))
	{
		points(pc.grid.return100W[i], return100LOGcW.glm$fitted.values[i], cex=0.5)
		points(pc.grid.return100W[i], (return100LOGcW.glm$fitted.values[i]+CI.z*pred100LOGW$se.fit[i]), cex = 0.4, col="red")
		points(pc.grid.return100W[i], (return100LOGcW.glm$fitted.values[i]-CI.z*pred100LOGW$se.fit[i]), cex = 0.4, col="red")

	}
dev.off() 


###########################################################################################
###########################################################################################
###########################################################################################


return100LOGelevC.glm <- glm(log(pc.return100) ~ pc.grid.return100 + pc.elev) 
pred100LOGElev <- predict(return100LOGelevC.glm, se.fit=TRUE)
out100LOGElev <- c(return100LOGelevC.glm$coefficients, return100LOGelevC.glm$aic, return100LOGelevC.glm$deviance,
length(return100LOGelevC.glm$residuals))
#setwd(wdMP)
write.table(t(out100LOGElev), file="PointVSGrid/WINTER/Models/GLM/100Ret/ModelParams/100RetDJFLOGPtElevModel.csv", col.names=param.headE)

#setwd(wdPNG)

png(filename="PointVSGrid/WINTER/Models/GLM/100Ret/PNG/100RetLogPtElevResidDJF.png")
plot(pc.grid.return100, return100LOGelevC.glm$residuals, main="100-Yr Ret LogPt+Elev Resid DJF", cex=0.4)
dev.off()

png(filename="PointVSGrid/WINTER/Models/GLM/100Ret/PNG/100RetLogPtElevStuResidDJF.png")
plot(pc.grid.return100, rstudent(return100LOGelevC.glm), main="100-Yr Ret LogPt+Elev StuResid DJF", cex=0.4)
dev.off()



png(filename="PointVSGrid/WINTER/Models/GLM/100Ret/PNG/100RetLogPtElevFitDJF.png")
plot(pc.grid.return100, log(pc.return100), main="100-Yr Ret LogPt+Elev Fitted Values DJF", cex=0.4, col="green")
for(i in 1:length(p.TEMP2[,1]))
	{
		points(pc.grid.return100[i], return100LOGelevC.glm$fitted.values[i], cex=0.5)
		points(pc.grid.return100[i], (return100LOGelevC.glm$fitted.values[i]+CI.z*pred100LOGElev$se.fit[i]), cex = 0.4, col="red")
		points(pc.grid.return100[i], (return100LOGelevC.glm$fitted.values[i]-CI.z*pred100LOGElev$se.fit[i]), cex = 0.4, col="red")

	} 
dev.off()

########################################################################
########################################################################
########################################################################
#Clean studentized residuals magnitude greater than 2
pc.return.weightE <- matrix(1, length(p.TEMP2[,1]), 1)

for(i in 1:length(p.TEMP2[,1]))  #This takes ~10 minutes to run
	{
		if(rstudent(return100LOGelevC.glm)[i] > 2 || rstudent(return100LOGelevC.glm)[i] < -2)
		{pc.return.weightE[i] <- 0}

	}

pc.grid.return100WE <- NULL
pc.return100WE <- NULL
pc.elevWE <- NULL

for(i in 1:length(p.TEMP2[,1]))  
	{
		if(pc.return.weightE[i] != 0)
		{pc.grid.return100WE <- rbind(pc.grid.return100WE, pc.grid.return100[i]);
		pc.return100WE <- rbind(pc.return100WE, pc.return100[i]);
		pc.elevWE <- rbind(pc.elevWE, pc.elev[i])}

	}
#setwd("/Users/shamseld/Desktop/NCAR/PointVSGrid/WINTER")
write.table(cbind(pc.grid.return100WE, pc.return100WE, pc.elevWE), file="PointVSGrid/WINTER/PDJF100WE.csv")


return100LOGelevCW.glm <- glm(log(pc.return100WE) ~ pc.grid.return100WE + pc.elevWE)  
pred100LOGElevW <- predict(return100LOGelevCW.glm, se.fit=TRUE)
out100LOGElevW <- c(return100LOGelevCW.glm$coefficients, return100LOGelevCW.glm$aic, return100LOGelevCW.glm$deviance,
length(return100LOGelevCW.glm$residuals))
#setwd(wdMP)
write.table(t(out100LOGElevW), file="PointVSGrid/WINTER/Models/GLM/100Ret/ModelParams/100RetDJFLOGPtElevModelW.csv", col.names=param.headE)

#setwd(wdPNG)

png(filename="PointVSGrid/WINTER/Models/GLM/100Ret/PNG/100RetLogPtElevResidDJFW.png")
plot(pc.grid.return100WE, return100LOGelevCW.glm$residuals, main="100-Yr Ret LogPt+Elev Resid DJFW", cex=0.4)
dev.off()

png(filename="PointVSGrid/WINTER/Models/GLM/100Ret/PNG/100RetLogPtElevStuResidDJFW.png")
plot(pc.grid.return100WE, rstudent(return100LOGelevCW.glm), main="100-Yr Ret LogPt+Elev StuResid DJFW", cex=0.4)
dev.off()


png(filename="PointVSGrid/WINTER/Models/GLM/100Ret/PNG/100RetLogPtElevFitDJFW.png")
plot(pc.grid.return100, log(pc.return100), main="100-Yr Ret LogPt+Elev Fitted Values DJF-W", cex=0.4, col="green")
for(i in 1:length(p.TEMP2[,1]))
	{
		points(pc.grid.return100[i], return100LOGelevCW.glm$fitted.values[i], cex=0.5)
		points(pc.grid.return100[i], (return100LOGelevCW.glm$fitted.values[i]+CI.z*pred100LOGElevW$se.fit[i]), cex = 0.4, col="red")
		points(pc.grid.return100[i], (return100LOGelevCW.glm$fitted.values[i]-CI.z*pred100LOGElevW$se.fit[i]), cex = 0.4, col="red")

	} 
dev.off()

########################################################################
########################################################################
########################################################################

pc.return50 <- matrix(0, length(p.TEMP2[,1]), 1)
pc.grid.return50 <- matrix(0, length(p.TEMP2[,1]), 1)
n.return50 <- 50

#Set return level based on xi=0, xi!=0
for(i in 1:length(p.TEMP2[,1]))
	{
		if(pc.xi[i] != 0) 
			{ pc.return50[i] <- (pc.mu[i] + exp(pc.logpsi[i])*((n.return50^pc.xi[i] - 1)/pc.xi[i])) }
		if(pc.grid.xi[i] != 0)
 			{ pc.grid.return50[i] <- (pc.grid.mu[i] + exp(pc.grid.logpsi[i])*((n.return50^pc.grid.xi[i] - 1)/pc.grid.xi[i])) }
		
		if(pc.xi[i] == 0)
			{ pc.return50[i] <- (pc.mu[i] + exp(pc.logpsi[i])*log(n.return50)) }
 		if(pc.grid.xi[i] == 0)
 			{ pc.grid.return50[i] <- (pc.grid.mu[i] + exp(pc.grid.logpsi[i])*log(n.return50)) }
 		
	}



return50c.glm <- glm(pc.return50 ~ pc.grid.return50)  
pred50 <- predict(return50c.glm, se.fit=TRUE)
out50 <- c(return50c.glm$coefficients, return50c.glm$aic, return50c.glm$deviance, length(return50c.glm$residuals))
#setwd(wdMP50)
write.table(t(out50), file="PointVSGrid/WINTER/Models/GLM/50Ret/ModelParams/50RetDJFModel.csv", col.names=param.head)

#setwd(wdPNG50)
png(filename="PointVSGrid/WINTER/Models/GLM/50Ret/PNG/50RetResidDJF.png", pointsize=6)
plot(pc.grid.return50, return50c.glm$residuals, main="50-Yr Ret Residuals - DJF", cex=0.5)
dev.off()


png(filename="PointVSGrid/WINTER/Models/GLM/50Ret/PNG/50RetStuResidDJF.png")
plot(pc.grid.return50, rstudent(return50c.glm), main="50-Yr Ret Student Residuals - DJF", cex=0.5)
dev.off()

png(filename="PointVSGrid/WINTER/Models/GLM/50Ret/PNG/50RetFitDJF.png")
plot(pc.grid.return50, pc.return50, main="50-Yr Ret Fitted Values", cex=0.5, col="green")
for(i in 1:length(p.TEMP2[,1]))
	{
		points(pc.grid.return50[i], return50c.glm$fitted.values[i], cex=0.4)
		points(pc.grid.return50[i], (return50c.glm$fitted.values[i]+CI.z*pred50$se.fit[i]), cex = 0.4, col="red")
		points(pc.grid.return50[i], (return50c.glm$fitted.values[i]-CI.z*pred50$se.fit[i]), cex = 0.4, col="red")
	}
dev.off() 



return50LOGc.glm <- glm(log(pc.return50) ~ pc.grid.return50) 
pred50LOG <- predict(return50LOGc.glm, se.fit=TRUE)
out50LOG <- c(return50LOGc.glm$coefficients, return50LOGc.glm$aic, return50LOGc.glm$deviance, length(return50LOGc.glm$residuals))
#setwd(wdMP50)
write.table(t(out50LOG), file="PointVSGrid/WINTER/Models/GLM/50Ret/ModelParams/50RetDJFLOGPtModel.csv", col.names=param.head)

#setwd(wdPNG50)
png(filename="PointVSGrid/WINTER/Models/GLM/50Ret/PNG/50RetLogPtResidDJF.png")
plot(pc.grid.return50, return50LOGc.glm$residuals, main="50-Yr Ret LogPt Resid - DJF", cex=0.4)
dev.off() 

png(filename="PointVSGrid/WINTER/Models/GLM/50Ret/PNG/50RetLogPtStuResidDJF.png")
plot(pc.grid.return50, rstudent(return50LOGc.glm), main="50-Yr Ret LogPt StuResid - DJF", cex=0.4)
dev.off() 

png(file="PointVSGrid/WINTER/Models/GLM/50Ret/PNG/50RetLogPtFitDJF.png")
plot(pc.grid.return50, log(pc.return50), main="50-Yr Return  LogPt Fitted Values - DJF", cex=0.4, col="green")
for(i in 1:length(p.TEMP2[,1]))
	{
		points(pc.grid.return50[i], return50LOGc.glm$fitted.values[i], cex=0.5)
		points(pc.grid.return50[i], (return50LOGc.glm$fitted.values[i]+CI.z*pred50LOG$se.fit[i]), cex = 0.4, col="red")
		points(pc.grid.return50[i], (return50LOGc.glm$fitted.values[i]-CI.z*pred50LOG$se.fit[i]), cex = 0.4, col="red")
	}
dev.off() 


#Clean studentized residuals magnitude greater than 2
#Weight based on studentized residuals
pc.return.weight <- matrix(1, length(p.TEMP2[,1]), 1)

for(i in 1:length(p.TEMP2[,1]))  #This takes ~10 minutes to run
	{
		if(rstudent(return50LOGc.glm)[i] > 2 || rstudent(return50LOGc.glm)[i] < -2)
		{pc.return.weight[i] <- 0}

	}

pc.grid.return50W <- NULL
pc.return50W <- NULL
for(i in 1:length(p.TEMP2[,1]))  
	{
		if(pc.return.weight[i] != 0)
		{pc.grid.return50W <- rbind(pc.grid.return50W, pc.grid.return50[i]);
		pc.return50W <- rbind(pc.return50W, pc.return50[i])}

	}

#setwd("/Users/shamseld/Desktop/NCAR/PointVSGrid/WINTER")
write.table(cbind(pc.grid.return50W, pc.return50W), file="PointVSGrid/WINTER/PDJF50W.csv")


return50LOGcW.glm <- glm(log(pc.return50W) ~ pc.grid.return50W) 
pred50LOGW <- predict(return50LOGcW.glm, se.fit=TRUE, )
out50LOGW <- c(return50LOGcW.glm$coefficients, return50LOGcW.glm$aic, return50LOGcW.glm$deviance,
length(return50LOGcW.glm$residuals))
#setwd(wdMP50)
write.table(t(out50LOGW), file="PointVSGrid/WINTER/Models/GLM/50Ret/ModelParams/50RetDJFLOGPtModelW.csv", col.names=param.head)

#setwd(wdPNG50)
png(file="PointVSGrid/WINTER/Models/GLM/50Ret/PNG/50RetLogPtStuResidDJFW.png")
plot(pc.grid.return50W, rstudent(return50LOGcW.glm), main="50-Yr Ret LogPt (W) Student Resid - DJF", cex=0.4)
dev.off() 

png(file="PointVSGrid/WINTER/Models/GLM/50Ret/PNG/50RetLogPtFitDJFW.png")
plot(pc.grid.return50W, log(pc.return50W), main="50-Yr Ret LogPt (W) Fitted Values - DJF", cex=0.4, col="green")
for(i in 1:length(p.TEMP2[,1]))
	{
		points(pc.grid.return50W[i], return50LOGcW.glm$fitted.values[i], cex=0.5)
		points(pc.grid.return50W[i], (return50LOGcW.glm$fitted.values[i]+CI.z*pred50LOGW$se.fit[i]), cex = 0.4, col="red")
		points(pc.grid.return50W[i], (return50LOGcW.glm$fitted.values[i]-CI.z*pred50LOGW$se.fit[i]), cex = 0.4, col="red")

	}
dev.off() 



return50LOGelevC.glm <- glm(log(pc.return50) ~ pc.grid.return50 + pc.elev) 
pred50LOGElev <- predict(return50LOGelevC.glm, se.fit=TRUE)
out50LOGElev <- c(return50LOGelevC.glm$coefficients, return50LOGelevC.glm$aic, return50LOGelevC.glm$deviance,
length(return50LOGelevC.glm$residuals))
#setwd(wdMP50)
write.table(t(out50LOGElev), file="PointVSGrid/WINTER/Models/GLM/50Ret/ModelParams/50RetDJFLOGPtElevModel.csv", col.names=param.headE)

#setwd(wdPNG50)

png(filename="PointVSGrid/WINTER/Models/GLM/50Ret/PNG/50RetLogPtElevResidDJF.png")
plot(pc.grid.return50, return50LOGelevC.glm$residuals, main="50-Yr Ret LogPt+Elev Resid DJF", cex=0.4)
dev.off()

png(filename="PointVSGrid/WINTER/Models/GLM/50Ret/PNG/50RetLogPtElevStuResidDJF.png")
plot(pc.grid.return50, rstudent(return50LOGelevC.glm), main="50-Yr Ret LogPt+Elev StuResid DJF", cex=0.4)
dev.off()



png(filename="PointVSGrid/WINTER/Models/GLM/50Ret/PNG/50RetLogPtElevFitDJF.png")
plot(pc.grid.return50, log(pc.return50), main="50-Yr Ret LogPt+Elev Fitted Values DJF", cex=0.4, col="green")
for(i in 1:length(p.TEMP2[,1]))
	{
		points(pc.grid.return50[i], return50LOGelevC.glm$fitted.values[i], cex=0.5)
		points(pc.grid.return50[i], (return50LOGelevC.glm$fitted.values[i]+CI.z*pred50LOGElev$se.fit[i]), cex = 0.4, col="red")
		points(pc.grid.return50[i], (return50LOGelevC.glm$fitted.values[i]-CI.z*pred50LOGElev$se.fit[i]), cex = 0.4, col="red")

	} 
dev.off()

########################################################################
########################################################################
########################################################################
#Clean studentized residuals magnitude greater than 2
pc.return.weightE <- matrix(1, length(p.TEMP2[,1]), 1)

for(i in 1:length(p.TEMP2[,1]))  #This takes ~10 minutes to run
	{
		if(rstudent(return50LOGelevC.glm)[i] > 2 || rstudent(return50LOGelevC.glm)[i] < -2)
		{pc.return.weightE[i] <- 0}

	}

pc.grid.return50WE <- NULL
pc.return50WE <- NULL
pc.elevWE <- NULL

for(i in 1:length(p.TEMP2[,1]))  
	{
		if(pc.return.weightE[i] != 0)
		{pc.grid.return50WE <- rbind(pc.grid.return50WE, pc.grid.return50[i]);
		pc.return50WE <- rbind(pc.return50WE, pc.return50[i]);
		pc.elevWE <- rbind(pc.elevWE, pc.elev[i])}

	}
#setwd("/Users/shamseld/Desktop/NCAR/PointVSGrid/WINTER")
write.table(cbind(pc.grid.return50WE, pc.return50WE, pc.elevWE), file="PointVSGrid/WINTER/PDJF50WE.csv")




return50LOGelevCW.glm <- glm(log(pc.return50WE) ~ pc.grid.return50WE + pc.elevWE)  
pred50LOGElevW <- predict(return50LOGelevCW.glm, se.fit=TRUE)
out50LOGElevW <- c(return50LOGelevCW.glm$coefficients, return50LOGelevCW.glm$aic, return50LOGelevCW.glm$deviance,
length(return50LOGelevCW.glm$residuals))
#setwd(wdMP50)
write.table(t(out50LOGElevW), file="PointVSGrid/WINTER/Models/GLM/50Ret/ModelParams/50RetDJFLOGPtElevModelW.csv", col.names=param.headE)

#setwd(wdPNG50)

png(filename="PointVSGrid/WINTER/Models/GLM/50Ret/PNG/50RetLogPtElevResidDJFW.png")
plot(pc.grid.return50WE, return50LOGelevCW.glm$residuals, main="50-Yr Ret LogPt+Elev Resid DJFW", cex=0.4)
dev.off()

png(filename="PointVSGrid/WINTER/Models/GLM/50Ret/PNG/50RetLogPtElevStuResidDJFW.png")
plot(pc.grid.return50WE, rstudent(return50LOGelevCW.glm), main="50-Yr Ret LogPt+Elev StuResid DJFW", cex=0.4)
dev.off()


png(filename="PointVSGrid/WINTER/Models/GLM/50Ret/PNG/50RetLogPtElevFitDJFW.png")
plot(pc.grid.return50, log(pc.return50), main="50-Yr Ret LogPt+Elev Fitted Values DJF-W", cex=0.4, col="green")
for(i in 1:length(p.TEMP2[,1]))
	{
		points(pc.grid.return50[i], return50LOGelevCW.glm$fitted.values[i], cex=0.5)
		points(pc.grid.return50[i], (return50LOGelevCW.glm$fitted.values[i]+CI.z*pred50LOGElevW$se.fit[i]), cex = 0.4, col="red")
		points(pc.grid.return50[i], (return50LOGelevCW.glm$fitted.values[i]-CI.z*pred50LOGElevW$se.fit[i]), cex = 0.4, col="red")

	} 
dev.off()

#
#
#
#
#
#
#
#
#
#http://tolstoy.newcastle.edu.au/R/help/05/05/4251.html
# f.mlm <- lm(cbind(y1,y2)~z1+z2)  


ret50ret100LOG.lm <- lm(cbind(log(pc.return50),log(pc.return100)) ~ cbind(pc.grid.return50, pc.grid.return100))

pred50100LOG <- predict(ret50ret100LOG.lm, se.fit=TRUE)
out50100LOG <- c(ret50ret100LOG.lm$coefficients, length(ret50ret100LOG.lm$residuals))
#setwd(wdMP50100)
write.table(t(out50100LOG), file="PointVSGrid/WINTER/Models/LM/50Ret100Ret/ModelParams/50100RetDJFLOGPtModel.csv", col.names=param.head50100)

#setwd(wdPNG50100)
png(filename="PointVSGrid/WINTER/Models/LM/50Ret100Ret/PNG/50-100RetLogPtResidDJF.png")
plot(ret50ret100LOG.lm$fitted.values, ret50ret100LOG.lm$residuals, main="50 & 100-Yr Ret LogPt Resid - DJF", cex=0.4)
dev.off() 

png(filename="PointVSGrid/WINTER/Models/LM/50Ret100Ret/PNG/50-100RetLogPtStuResidDJF.png")
plot(ret50ret100LOG.lm$fitted.values, rstudent(ret50ret100LOG.lm), main="50 & 100-Yr Ret LogPt StuResid - DJF", cex=0.4)
dev.off() 

png(file="PointVSGrid/WINTER/Models/LM/50Ret100Ret/PNG/50Ret100RetLogPtFitDJF.png")
plot(pc.grid.return100, log(pc.return100), main="50& 100-Yr Return  LogPt Fitted Values - DJF", cex=0.3, col="green")
for(i in 1:length(pc.grid.return50))
	{
		points(pc.grid.return50[i], log(pc.return50[i]), cex=0.3, col="blue")
		points(pc.grid.return50[i], ret50ret100LOG.lm$fitted.values[i,1], cex=0.8)
		points(pc.grid.return100[i], ret50ret100LOG.lm$fitted.values[i,2], cex=0.8)

	}

dev.off() 

ret50ret100LOGE.lm <- lm(cbind(log(pc.return50),log(pc.return100)) ~ cbind(pc.grid.return50, pc.grid.return100)+pc.elev)

pred50100LOGE <- predict(ret50ret100LOGE.lm, se.fit=TRUE)
out50100LOGE <- c(ret50ret100LOGE.lm$coefficients, length(ret50ret100LOGE.lm$residuals))
#setwd(wdMP50100)
write.table(t(out50100LOGE), file="PointVSGrid/WINTER/Models/LM/50Ret100Ret/ModelParams/50100RetDJFLOGPtElevModel.csv", col.names=param.head50100E1)

#setwd(wdPNG50100)
png(filename="PointVSGrid/WINTER/Models/LM/50Ret100Ret/PNG/50-100RetLogPtElevResidDJF.png")
plot(ret50ret100LOGE.lm$fitted.values, ret50ret100LOGE.lm$residuals, main="50 & 100-Yr Ret LogPt+Elev Resid - DJF", cex=0.4)
dev.off() 

png(filename="PointVSGrid/WINTER/Models/LM/50Ret100Ret/PNG/50-100RetLogPtElevStuResidDJF.png")
plot(ret50ret100LOGE.lm$fitted.values, rstudent(ret50ret100LOGE.lm), main="50 & 100-Yr Ret LogPt+Elev StuResid - DJF", cex=0.4)
dev.off() 

png(file="PointVSGrid/WINTER/Models/LM/50Ret100Ret/PNG/50Ret100RetLogPtElevFitDJF.png")
plot(pc.grid.return100, log(pc.return100), main="50& 100-Yr Return  LogPt+Elev Fitted Values - DJF", cex=0.3, col="green")
for(i in 1:length(pc.grid.return50))
	{
		points(pc.grid.return50[i], log(pc.return50[i]), cex=0.3, col="blue")
		points(pc.grid.return50[i], ret50ret100LOGE.lm$fitted.values[i,1], cex=0.8)
		points(pc.grid.return100[i], ret50ret100LOGE.lm$fitted.values[i,2], cex=0.8)

	}

dev.off() 

#Clean studentized residuals magnitude greater than 2
#Weight based on studentized residuals
Sys.time()
pc.return.weight50100 <- matrix(1, length(ret50ret100LOG.lm$fitted.values), 1)

for(i in 1:length(ret50ret100LOG.lm$fitted.values))  #This takes ~18 minutes to run
	{
		if(rstudent(ret50ret100LOG.lm)[i] > 2 || rstudent(ret50ret100LOG.lm)[i] < -2)
		{pc.return.weight50100[i] <- 0}

	}
pc.grid.return50100W <- NULL
pc.return50100W <- NULL
for(i in 1:(length(pc.return.weight50100)/2))
	{
		if(pc.return.weight50100[i] != 0)
		{pc.grid.return50100W <- rbind(pc.grid.return50100W, pc.grid.return50[i]);
		pc.return50100W <- rbind(pc.return50100W, pc.return50[i])}
	}

for(i in (length(pc.return.weight50100)/2 + 1):length(pc.return.weight50100))
	{
		if(pc.return.weight50100[i] != 0)
		{pc.grid.return50100W <- rbind(pc.grid.return50100W, pc.grid.return100[(i-length(pc.return.weight50100)/2)]);
		pc.return50100W <- rbind(pc.return50100W, pc.return100[(i-length(pc.return.weight50100)/2)])}
	}

Sys.time()


#setwd("/Users/shamseld/Desktop/NCAR/PointVSGrid/WINTER")
write.table(cbind(pc.grid.return50100W, pc.return50100W), file="PointVSGrid/WINTER/PDJF50100W.csv")

ret50ret100LOGW.lm <- lm(cbind(log(pc.return50100W[1:(length(pc.return50100W)/2)]),
	log(pc.return50100W[(length(pc.return50100W)/2 + 1):length(pc.return50100W)])) ~ 
	cbind(pc.grid.return50100W[1:(length(pc.return50100W)/2)], 	pc.grid.return50100W[(length(pc.return50100W)/2 + 1):length(pc.return50100W)]))

pred50100LOGW <- predict(ret50ret100LOGW.lm, se.fit=TRUE)
out50100LOGW <- c(ret50ret100LOGW.lm$coefficients, length(ret50ret100LOGW.lm$fitted.values))
#setwd(wdMP50100)
write.table(t(out50100LOGW), file="PointVSGrid/WINTER/Models/LM/50Ret100Ret/ModelParams/50100RetDJFLOGPtModelW.csv", col.names=param.head50100)

#setwd(wdPNG50100)
png(filename="PointVSGrid/WINTER/Models/LM/50Ret100Ret/PNG/50-100RetLogPtResidDJFW.png")
plot(ret50ret100LOGW.lm$fitted.values, ret50ret100LOGW.lm$residuals, main="50 & 100-Yr Ret LogPt (W) Resid - DJF", cex=0.4)
dev.off() 

png(filename="PointVSGrid/WINTER/Models/LM/50Ret100Ret/PNG/50-100RetLogPtStuResidDJFW.png")
plot(ret50ret100LOGW.lm$fitted.values, rstudent(ret50ret100LOGW.lm), main="50 & 100-Yr Ret LogPt (W) StuResid - DJF", cex=0.4)
dev.off() 

png(file="PointVSGrid/WINTER/Models/LM/50Ret100Ret/PNG/50Ret100RetLogPtFitDJFW.png")
plot(pc.grid.return50100W, log(pc.return50100W), main="50& 100-Yr Return  LogPt (W) Fitted - DJF", cex=0.3, col="green")
for(i in 1:length(pc.return.weight50100))
	{
		points(pc.grid.return50100W[i], ret50ret100LOGW.lm$fitted.values[i], cex=0.8)

	}

dev.off() 


#Clean studentized residuals magnitude greater than 2
#Weight based on studentized residuals
Sys.time()
pc.return.weight50100E <- matrix(1, length(ret50ret100LOG.lm$fitted.values), 1)

for(i in 1:length(ret50ret100LOGE.lm$fitted.values))  #This takes ~8 minutes to run
	{
		if(rstudent(ret50ret100LOGE.lm)[i] > 2 || rstudent(ret50ret100LOGE.lm)[i] < -2)
		{pc.return.weight50100E[i] <- 0}

	}
pc.grid.return50100WE <- NULL
pc.return50100WE <- NULL
pc.elev50100WE <- NULL

for(i in 1:(length(pc.return.weight50100E)/2))
	{
		if(pc.return.weight50100E[i] != 0)
		{pc.grid.return50100WE <- rbind(pc.grid.return50100WE, pc.grid.return50[i]);
		pc.return50100WE <- rbind(pc.return50100WE, pc.return50[i]);
		pc.elev50100WE <- rbind(pc.elev50100WE, pc.elev[i])}
	}

for(i in (length(pc.return.weight50100)/2 + 1):length(pc.return.weight50100))
	{
		if(pc.return.weight50100E[i] != 0)
		{pc.grid.return50100WE <- rbind(pc.grid.return50100WE, pc.grid.return100[(i-length(pc.return.weight50100)/2)]);
		pc.return50100WE <- rbind(pc.return50100WE, pc.return100[(i-length(pc.return.weight50100)/2)]);
		pc.elev50100WE <- rbind(pc.elev50100WE, pc.elev[(i-length(pc.return.weight50100)/2) ] ) }
	}

Sys.time()


#setwd("/Users/shamseld/Desktop/NCAR/PointVSGrid/WINTER")
write.table(cbind(pc.grid.return50100WE, pc.return50100WE, pc.elev50100WE), file="PointVSGrid/WINTER/PDJF50100WE.csv")


ret50ret100LOGWE.lm <- lm(cbind(log(pc.return50100WE[1:(length(pc.return50100WE)/2)]),
	log(pc.return50100WE[(length(pc.return50100WE)/2 + 1):length(pc.return50100WE)])) ~ 
	cbind(pc.grid.return50100WE[1:(length(pc.return50100WE)/2)], 	pc.grid.return50100WE[(length(pc.return50100WE)/2 + 1):length(pc.return50100WE)]) 
	+ cbind(pc.elev50100WE[1:(length(pc.return50100WE)/2)], pc.elev50100WE[(length(pc.return50100WE)/2 + 1):length(pc.return50100WE)]))

pred50100LOGWE <- predict(ret50ret100LOGWE.lm, se.fit=TRUE)
out50100LOGWE <- c(ret50ret100LOGWE.lm$coefficients, length(ret50ret100LOGWE.lm$fitted.values))
#setwd(wdMP50100)
write.table(t(out50100LOGWE), file="PointVSGrid/WINTER/Models/LM/50Ret100Ret/ModelParams/50100RetDJFLOGPtElevModelW.csv",
col.names=param.head50100E2)

#setwd(wdPNG50100)
png(filename="PointVSGrid/WINTER/Models/LM/50Ret100Ret/PNG/50-100RetLogPtElevResidDJFW.png")
plot(ret50ret100LOGWE.lm$fitted.values, ret50ret100LOGWE.lm$residuals, main="50 & 100-Yr Ret LogPt+Elev (W) Resid - DJF", cex=0.4)
dev.off() 

png(filename="PointVSGrid/WINTER/Models/LM/50Ret100Ret/PNG/50-100RetLogPtElevStuResidDJFW.png")
plot(ret50ret100LOGWE.lm$fitted.values, rstudent(ret50ret100LOGWE.lm), main="50 & 100-Yr Ret LogPt+Elev (W) StuResid - DJF", cex=0.4)
dev.off() 

png(file="PointVSGrid/WINTER/Models/LM/50Ret100Ret/PNG/50Ret100RetLogPtElevFitDJFW.png")
plot(pc.grid.return50100WE, log(pc.return50100WE), main="50& 100-Yr Return  LogPt+Elev (W) Fitted - DJF", cex=0.3, col="green")
for(i in 1:length(pc.return50100WE))
	{
		points(pc.grid.return50100WE[i], ret50ret100LOGWE.lm$fitted.values[i], cex=0.8)

	}

dev.off() 

EndTime <- Sys.time()

StartTime
EndTime

###########################################################################################
###########################################################################################
###########################################################################################