######################################################################################################
#
#		R CODE FOR ANALYSIS OF US PRECIP DATA BASED ON PRECIP YEAR OCT - SEP
#		USES ALL 20 CITIES OR SUBSET OF LONGEST RECORDS
#
######################################################################################################

rm(list=ls(all=TRUE))			# remove all objects in memory
options(scipen=999)			# forbid printing in scientific notation
library(readxl)	
library(lattice)
library(fracdiff)
library(maps)


######################################################################################################
###	SET PARAMETERS
startyear = 1978						# Choose starting year (earliest 1889 or 1872)
endyear = 2018							# Choose end year (note precip year starts Oct before)
seasonal =0							# Annual (0) or seasonal (1) analysis
season.select = 4						# Season selected winter = 1, spring = 2, summer = 3, fall = 4
fewlong = 0							# Do analysis on smaller full-length set
pentad.graph = 0						# Do 5-year event graph
vcov = 1							# Compute VF matrix and write to file (alternative is to read from file)
fd = 0								# Compute fractional difference parms
######################################################################################################




######################################################################################################
###	FUNCTION OMEGA1 takes matrix x (residuals) with number of series K>1 and length N
###	then generates VF2005 omega matrix (no breaks)

omega1 = function(x)	{
	K = ncol(x)
	N = nrow(x)
	d = c(1, rep( 0 , times=N-1))				# create N-length vector of 1 followed by N-1 zeroes
	x = t(x)						# transpose residual matrix
	S= x[1:K,1] 						# S = first col of u 
	OMEGA1 = S %*% t(S)					# First accum entry for Omega matrix
	for (i in 2:N)		{				# Loop to accumulate other Omega=SS' terms
		d[i]=1						# Vector to add 1st i cols of u (sets next element to 1)
		S = x %*% d					# forms Nx1 vector with sum of 1st i cols of u
		OMEGA1 = OMEGA1 + S %*% t(S)			# accumulates partial sums
				}
	OMEGA1 = 2*OMEGA1 / (N^2)				# Last step
	return(OMEGA1)
				}

###	FUNCTION OMEGA0 takes column x (residuals) with number of series K=1 and length N
###	then generates VF2005 omega matrix (no breaks)
omega0 = function(x)		{
	K = 1
	N = length(r_t)
	d = c(1, rep( 0 , times=N-1))					# create T-length vector of 1 followed by N-1 zeroes
	x = t(x)							# transpose residual matrix
	S= x[1:K,1] 							# S = first col of u 
	OMEGA1 = S %*% t(S)						# First accum entry for Omega matrix
	for (i in 2:N)		{					# Loop to accumulate other Omega=SS' terms
		d[i]=1							# Vector to add 1st i cols of u (sets next element to 1)
		S = x %*% d						# forms Nx1 vector with sum of 1st i cols of u
		OMEGA1 = OMEGA1 + S %*% t(S)				# accumulates partial sums
				}
	OMEGA1 = 2*OMEGA1 / (N^2)					# Last step
	return(OMEGA1)
				}
######################################################################################################





######################################################################################################
###	GET LOCATION NAMES & MAKE SPATIAL WEIGHTING VECTORS
names = read.csv("data/names.csv", head=FALSE, sep=",")		# Location names; switch to few long if indicated
	S_pc = as.matrix(c(1,1,1,1,1,1,1,1,1,1,0,0,0,0,0,0,0,0,0,0)/10)
	S_se = as.matrix(c(0,0,0,0,0,0,0,0,0,0,1,1,1,1,1,1,1,1,1,1)/10)
	S_all= as.matrix(c(1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1)/20)

if (fewlong==1) {
	names = read.csv("data/names2.csv", head=FALSE, sep=",")
		S_pc = as.matrix(c(1,1,1,1,1,0,0,0,0,0,0)/5)
		S_se = as.matrix(c(0,0,0,0,0,1,1,1,1,1,1)/6)
		S_all= as.matrix(c(1,1,1,1,1,1,1,1,1,1,1)/11)
		}

loc_code = names$V1						# Number of names and location codes
loc_name = names$V2
region = names$V3
lat = names$V4
lon = names$V5
loc_n	 = length(loc_code)

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




######################################################################################################
###	CONSTRUCT DATA SETS 

downpour = 25.4							# Cut-off (in mm) for daily downpour

#	MATRICES TO STORE AND PROCESS RESULTS 
nyears1 = endyear - startyear +1				# Number of years in drawn sample
nyears = nyears1-1						# Number of years in used sample
year = seq(startyear, (endyear-1))				# Sequence of years in used sample
tt = c(1:nyears)						# Simple time trend 1,2,...,T

#	DAY AND SEASON INDICATOR
day = c(1:366)							# day index
season = rep(1, times = 366)					# seasons: fall = 1, winter = 2, spring = 3, summer = 4
season[day >= 93 & day < 184] = 2
season[day >= 184 & day < 275] = 3
season[day >= 275] = 4

ndays = 366							# Number of days kept in year
if (seasonal == 1 & season.select == 1) {ndays = 92}		# 
if (seasonal == 1 & season.select == 2) {ndays = 91}
if (seasonal == 1 & season.select == 3) {ndays = 91}
if (seasonal == 1 & season.select == 4) {ndays = 92}

daily = matrix(0, ncol = loc_n, nrow = ndays*nyears)		# (BIG) matrix of daily values
daily_2day = matrix(0, ncol = loc_n, nrow = ndays*nyears)	# (BIG) matrix of 2-day totals

ann3_avg = matrix(0, ncol = loc_n, nrow = nyears, dimnames=list(year,loc_name) )	# matrix of annual values
ann4_med = matrix(0, ncol = loc_n, nrow = nyears, dimnames=list(year,loc_name) )
ann5_var = matrix(0, ncol = loc_n, nrow = nyears, dimnames=list(year,loc_name) )
ann6_max = matrix(0, ncol = loc_n, nrow = nyears, dimnames=list(year,loc_name) )
ann7_dpr = matrix(0, ncol = loc_n, nrow = nyears, dimnames=list(year,loc_name) )
ann_1d95 = matrix(0, ncol = loc_n, nrow = nyears, dimnames=list(year,loc_name) )
ann_1d99 = matrix(0, ncol = loc_n, nrow = nyears, dimnames=list(year,loc_name) )
ann_2d95 = matrix(0, ncol = loc_n, nrow = nyears, dimnames=list(year,loc_name) )
ann_2d99 = matrix(0, ncol = loc_n, nrow = nyears, dimnames=list(year,loc_name) )

tabcols = c("Start Year", "Mean", "Median", "Variance", "Max", "Trend", "CI-lower", "CI-upper", "VF Score")
Table1_1day_total = matrix(0, ncol=9, nrow=loc_n, dimnames=list(loc_name,tabcols) )
Table2_2day_total = matrix(0, ncol=9, nrow=loc_n, dimnames=list(loc_name,tabcols) )

tabcols = c("Trend", "CI-lower", "CI-upper", "VF Score")	# Results Tables
tabrows = c(paste(loc_name), "Pacific", "South East", "All")
Table3_avg = matrix(0, ncol=4, nrow=loc_n+3, dimnames=list(tabrows,tabcols) )
Table4_med = matrix(0, ncol=4, nrow=loc_n+3, dimnames=list(tabrows,tabcols) )
Table5_var = matrix(0, ncol=4, nrow=loc_n+3, dimnames=list(tabrows,tabcols) )
Table6_max = matrix(0, ncol=4, nrow=loc_n+3, dimnames=list(tabrows,tabcols) )
Table6_max2 = matrix(0, ncol=4, nrow=loc_n+3, dimnames=list(tabrows,tabcols) )
Table7_dpr = matrix(0, ncol=4, nrow=loc_n+3, dimnames=list(tabrows,tabcols) )
Table8_1day_99 = matrix(0, ncol=4, nrow=loc_n+3, dimnames=list(tabrows,tabcols) )
Table9_2day_99 = matrix(0, ncol=4, nrow=loc_n+3, dimnames=list(tabrows,tabcols) )

Reg_avg = c(" 0 "," 0 "," 0 ")					# Regional Result codes
Reg_med = c(" 0 "," 0 "," 0 ")	
Reg_var = c(" 0 "," 0 "," 0 ")	
Reg_max = c(" 0 "," 0 "," 0 ")	
Reg_dpr = c(" 0 "," 0 "," 0 ")	
Reg_1dx = c(" 0 "," 0 "," 0 ")	
Reg_2dx = c(" 0 "," 0 "," 0 ")	
Region = data.frame(Reg_avg,Reg_med,Reg_var,Reg_max,Reg_dpr,Reg_1dx,Reg_2dx, row.names = c("Pacific", "SouthEast", "Total"), 
		stringsAsFactors=FALSE)

p5 =matrix(0,nrow= nyears/5, ncol=loc_n)				# counts of pentad events
p50=matrix(0,nrow= 75, ncol=loc_n)







#	BUILD MAIN DATA MATRICES
#	START LOOP OVER LOCATIONS
for (j in 1:loc_n) 	{

#	EXTRACT DATA FROM EXCEL SPREADSHEET = dat1
firstyear = read_excel("data/dly_20190208_v_e.xlsx", sheet=paste(loc_code[j]), range = "D5", col_names=FALSE)
cols.to.read = endyear-firstyear+1				# Read from 1st year up to endyear
dat1 = read_excel("data/dly_20190208_v_e.xlsx", sheet=paste(loc_code[j]), 
	range = anchored("D5", dim = c(367, cols.to.read)), col_names = FALSE)
dat1 = data.frame(dat1)

#	SET NON-LEAP YEARS TO NA
for (i in 1:as.numeric(cols.to.read)) { dat1[61,i][dat1[61,i]==-99] = NA }

#	SELECT COLUMNS AFTER STARTYEAR = dat2
dat2 = subset(dat1, select = (dat1[1,]>= startyear))
rm(dat1)

#	RECTANGULAR DATA BY YEAR (LOCATION j) BUT WITHOUT YEAR HEADING = dat3
if (ncol(dat2) != nyears1) {cat("INCORRECT NUMBER OF YEARS READ FOR", loc_name[j], "\n")} 	# Check correct num years read
dat3 = dat2[2:367,]						# take daily data rows
dat3 = dat3/10							# convert to mm/day
rm(dat2)
n_miss = sum(dat3 < 0)

#	CONVERT TO PRECIP YEAR OCT - SEP = dat4, dat5, dat6, dpa_y
dat4 = stack(dat3, na.rm=TRUE)					# stack whole matrix into 1 long vector
big_N = length(dat4[,1])					# length of stacked vector
StartOct1 = 275							# trim up to 1st Sep 30 and from last Oct 1 on 
EndSep30 = big_N-92
dat5 = dat4[1:EndSep30,]					# remove last 92 obs (Oct 1 - Dec 31)
dat6 = dat5[StartOct1:EndSep30,]				# remove first 273 obs (Jan 1 - Sep 30)
dpa_y = matrix(dat6[,1], nrow=366, ncol=nyears)			# put into matrix dpa_y
rm(dat3, dat4, dat5, dat6)

#	USE SEASONAL RATHER THAN FULL-YEAR DATA IF SELECTED
dpa_ys = cbind(season[1:366], dpa_y)
if (seasonal == 1) {dpa_y = dpa_ys[season == season.select,2:(nyears+1)] }
rm(dpa_ys)

# 	ONE LONG SERIES "dpa" in matrix "daily"
dpa = stack(as.data.frame(dpa_y), na.rm=TRUE)[,1]		# make one long vector
big_N = length(dpa)						# number of obs
big_T = c(1:big_N)						# simple time trend
daily[,j] = dpa							# assign to matrix

# 	99th PERCENTILE CUT-OFF USING NON-ZERO DAYS
dpa2 = dpa
nz_dpa1 = dpa[dpa2>0]
cutoff99 = quantile(nz_dpa1, probs=0.99, na.rm=TRUE)		# 99th pctile cut-off
rm(dpa2)

#	2-DAY CUMULATIVES, PERCENTILE CUT-OFF & PENTAD EVENTS
lag_dpa = c(0,dpa[1:(big_N-1)])					# lag of dpa
dpa_2day = dpa+lag_dpa						# today's dpa + yesterday's
dpa_2day_y = matrix(dpa_2day, nrow=ndays, ncol=nyears)		# dpa_y made up of 2-day counts
dpa2 = dpa_2day							# compute 2-day percentiles
daily_2day[,j] = dpa_2day					# assign to matrix
nz_dpa2 = dpa_2day[dpa2>0]					#
cutoff99_2day = quantile(nz_dpa2, probs=0.99, na.rm=TRUE)	# 99th pctile cut-off
o_2day = order(dpa2, decreasing=TRUE)				# rank 2-day events
top_25 = dpa2[o_2day][25]					# pick 25th top
top_75 = dpa2[o_2day][75]					# pick 75th top
if(pentad.graph==1) 	{d2 = matrix(dpa2,nrow=5*366)		# organize daily 2 day data into pentads
	for (k in 1:ncol(d2))	{ p5[k,j] = sum(d2[,k] >= top_25, na.rm=TRUE)}	 			# count 5y events per pentad
			}
if(pentad.graph==1)	{
	for (k in 51:125) 	{p50[(k-50),j] = sum(dpa_2day_y[,(k-50):k] >=top_75, na.rm=TRUE)} 	# count 50y events per overlapping 50y period
			}

#	EXTRACT ANNUAL MEAN, VARIANCE AND MAX VALUES
ann3_avg[,j] = apply(dpa_y,2,mean, na.rm=TRUE)
ann5_var[,j] = apply(dpa_y,2,var, na.rm=TRUE)
ann6_max[,j] = apply(dpa_y,2,max, na.rm=TRUE)


#	ANNUAL NON-ZERO MEDIANS
for (k in 1:nyears) 	{ann4_med[k,j] = median(dpa_y[,k][dpa_y[,k]>0], na.rm=TRUE) }

#	ONE AND TWO-DAY ANNUAL EXCEEDANCE AND DOWNPOUR COUNTS
p_99 = c(0, times=nyears)
p_99_2d  = c(0, times=nyears)

for (i in 1:nyears) 	{
	ann_1d99[i,j] = sum(dpa_y[,i] >= cutoff99, na.rm=TRUE) 
	ann_2d99[i,j] = sum(dpa_2day_y[,i] >= cutoff99_2day, na.rm=TRUE) 
	 ann7_dpr[i,j] = sum(dpa_y[,i]>downpour, na.rm=TRUE) 
			}

#	PARTIALLY POPULATE TABLE 1
Table1_1day_total[j,1] = startyear
Table1_1day_total[j,2] = round(mean(dpa, na.rm=TRUE),2)
Table1_1day_total[j,3] = round(median(nz_dpa1, na.rm=TRUE),2)
Table1_1day_total[j,4] = round(var(dpa, na.rm=TRUE),2)
Table1_1day_total[j,5] = round(max(dpa, na.rm=TRUE),2)

#	PARTIALLY POPULATE TABLE 2
Table2_2day_total[j,1] = startyear
Table2_2day_total[j,2] = round(mean(dpa_2day, na.rm=TRUE),2)
Table2_2day_total[j,3] = round(median(nz_dpa2, na.rm=TRUE),2)
Table2_2day_total[j,4] = round(var(dpa_2day, na.rm=TRUE),2)
Table2_2day_total[j,5] = round(max(dpa_2day, na.rm=TRUE),2)

#	END DATA CONSTRUCTION LOOP
		}

cat("DATA SETS COMPLETED; STARTYEAR IS",startyear, "\n")

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



######################################################################################################
###	TREND ANALYSIS TO POPULATE TABLES

# 	TRENDS IN 1 & 2 DAY DAILY RECORDS - TABLES 1&2
K = loc_n
beta1 = lm(big_T ~ 1)						# partial out the t-terms
	ttilde1 = as.matrix(beta1$residuals)
	eta1 = sum( ttilde1^2 )					# partialling out term (no break)
reg_t1  = lm(daily ~ big_T)					# regress 1 day precip series on time trend T
	b_t1 = reg_t1$coefficients[2,]				# get slope coeffs
	r_t1 = reg_t1$residuals[,]				# Get residuals  
if (vcov == 1) {						# when Omega matrix is to be computed
	OMEGA1d=omega1(r_t1)
	if(fewlong==1) 	{write(OMEGA1d, file = "logfiles/omega1df.csv", ncol=loc_n, sep=",")} else
			{write(OMEGA1d, file = "logfiles/omega1d.csv", ncol=loc_n, sep=",")}
		}
	if(fewlong==1) 	{OMEGA1d=as.matrix(read.csv(file="logfiles/omega1df.csv", sep=",", head=FALSE))} else
			{OMEGA1d=as.matrix(read.csv(file="logfiles/omega1d.csv", sep=",", head=FALSE))}

reg_t2  = lm(daily_2day ~ big_T)				# regress 2 day precip series on time trend T
	b_t2 = reg_t2$coefficients[2,]				# get slope coeffs
	r_t2 = reg_t2$residuals[,]				# Get residuals  
if (vcov == 1) {						# when Omega matrix is to be computed
	OMEGA2d=omega1(r_t2)
	if(fewlong==1) 	{write(OMEGA2d, file = "logfiles/omega2df.csv", ncol=loc_n, sep=",")} else
			{write(OMEGA2d, file = "logfiles/omega2d.csv", ncol=loc_n, sep=",")}
		}
	if(fewlong==1) 	{OMEGA2d=as.matrix(read.csv(file="logfiles/omega2df.csv", sep=",", head=FALSE))} else
			{OMEGA2d=as.matrix(read.csv(file="logfiles/omega2d.csv", sep=",", head=FALSE))}


#	LOOP TO PUT TREND INFO IN TABLES 1 & 2 IN MM/YEAR (365 DAYS)
R=matrix(0, nrow=1, ncol=K)
for (i in 1:loc_n) 	{
	R[i] = 1
	VF = eta1 * (R %*% b_t1) %*% (solve( R %*% OMEGA1d %*% t(R) ) ) %*% (R %*% b_t1)
	Table1_1day_total[i,6] = round( 365*b_t1[i], 4)
	Table1_1day_total[i,7] = round( 365*(b_t1[i] - 6.482*abs( b_t1[i]/chol(VF))), 4)
	Table1_1day_total[i,8] = round( 365*(b_t1[i] + 6.482*abs( b_t1[i]/chol(VF))), 4)
	Table1_1day_total[i,9] = round(VF, 2)

	VF = eta1 * (R %*% b_t2) %*% (solve( R %*% OMEGA2d %*% t(R) ) ) %*% (R %*% b_t2)
	Table2_2day_total[i,6] = round( 365*b_t2[i], 4)
	Table2_2day_total[i,7] = round( 365*(b_t2[i] - 6.482*abs( b_t2[i]/chol(VF))), 4)
	Table2_2day_total[i,8] = round( 365*(b_t2[i] + 6.482*abs( b_t2[i]/chol(VF))), 4)
	Table2_2day_total[i,9] = round(VF, 2)

	R[i] = 0
			}
cat("TABLES 1 & 2 COMPLETED", "\n")

# 	TRENDS IN ANNUAL STATS
beta1 = lm(tt ~ 1)						# partial out the t-terms (tt is annual trend)
ttilde1 = as.matrix(beta1$residuals)
eta1 = sum( ttilde1^2 )						# partialling out term (no break)

#	AVG
#ann3_avg = ann3_avg						
reg_t  = lm(ann3_avg ~ tt)					# regress annual stat (all locations) on time trend
b_t = reg_t$coefficients[2,]					# get slope coeffs
r_t = reg_t$residuals[,]					# Get residuals  
OMEGA1=omega1(r_t)
R=matrix(0, nrow=1, ncol=loc_n)
for (j in 1:loc_n)	{					# j-loop for locations
	R[j] = 1
	VF = eta1 * (R %*% b_t) %*% (solve( R %*% OMEGA1 %*% t(R) ) ) %*% (R %*% b_t)
	Table3_avg[j,1] = round(b_t[j],4)
	Table3_avg[j,2] = round( (b_t[j] - 6.482*abs( b_t[j]/chol(VF))), 4)
	Table3_avg[j,3] = round( (b_t[j] + 6.482*abs( b_t[j]/chol(VF))), 4)
	Table3_avg[j,4] = round(VF,4)
	R[j] = 0
			}
# Pacific Region
	pc = ann3_avg %*% S_pc	
	pc_avg=pc 						# will need this later
	reg_t  = lm(pc ~ tt)
	b_t = reg_t$coefficients[2]
	r_t = reg_t$residuals
	OMEGA0=omega0(r_t)
	VF = eta1 * b_t^2 / OMEGA0
	Table3_avg[(loc_n+1),1] = round(b_t,4)
	Table3_avg[(loc_n+1),2] = round( (b_t - 6.482*abs(b_t/chol(VF))), 4)
	Table3_avg[(loc_n+1),3] = round( (b_t + 6.482*abs(b_t/chol(VF))), 4)
	Table3_avg[(loc_n+1),4] = round(VF, 4)
	if (b_t > 0 & VF > 41.53) {Region[1,1] = " + "}
	if (b_t < 0 & VF < 41.53) {Region[1,1] = "(0)"}
	if (b_t < 0 & VF > 41.53) {Region[1,1] = " - "}
# SouthEast Region
	se = ann3_avg %*% S_se
	se_avg=se 						# will need this later
	reg_t  = lm(se ~ tt)
	b_t = reg_t$coefficients[2]
	r_t = reg_t$residuals
	OMEGA0=omega0(r_t)
	VF = eta1 * b_t^2 / OMEGA0
	Table3_avg[(loc_n+2),1] = round(b_t,4)
	Table3_avg[(loc_n+2),2] = round( (b_t - 6.482*abs(b_t/chol(VF))), 4)
	Table3_avg[(loc_n+2),3] = round( (b_t + 6.482*abs(b_t/chol(VF))), 4)
	Table3_avg[(loc_n+2),4] = round(VF, 4)
	if (b_t > 0 & VF > 41.53) {Region[2,1] = " + "}
	if (b_t < 0 & VF < 41.53) {Region[2,1] = "(0)"}
	if (b_t < 0 & VF > 41.53) {Region[2,1] = " - "}
# All 
	all = ann3_avg %*% S_all
	all_avg = all						# will need this later
	reg_t  = lm(all ~ tt)
	b_t = reg_t$coefficients[2]
	r_t = reg_t$residuals
	OMEGA0=omega0(r_t)
	VF = eta1 * b_t^2 / OMEGA0
	Table3_avg[(loc_n+3),1] = round(b_t,4)
	Table3_avg[(loc_n+3),2] = round( (b_t - 6.482*abs(b_t/chol(VF))), 4)
	Table3_avg[(loc_n+3),3] = round( (b_t + 6.482*abs(b_t/chol(VF))), 4)
	Table3_avg[(loc_n+3),4] = round(VF, 4)
	if (b_t > 0 & VF > 41.53) {Region[3,1] = " + "}
	if (b_t < 0 & VF < 41.53) {Region[3,1] = "(0)"}
	if (b_t < 0 & VF > 41.53) {Region[3,1] = " - "}
#
cat("TABLE 3 (AVG) COMPLETE", "\n")

#	MED
reg_t  = lm(ann4_med ~ tt)					# regress annual stat (all locations) on time trend
b_t = reg_t$coefficients[2,]					# get slope coeffs
r_t = reg_t$residuals[,]					# Get residuals  
OMEGA1=omega1(r_t)
R=matrix(0, nrow=1, ncol=loc_n)
for (j in 1:loc_n)	{					# j-loop for locations
	R[j] = 1
	VF = eta1 * (R %*% b_t) %*% (solve( R %*% OMEGA1 %*% t(R) ) ) %*% (R %*% b_t)
	Table4_med[j,1] = round(b_t[j], 4)
	Table4_med[j,2] = round(b_t[j] - 6.482*abs( b_t[j]/chol(VF)), 4)
	Table4_med[j,3] = round(b_t[j] + 6.482*abs( b_t[j]/chol(VF)), 4)
	Table4_med[j,4] = round(VF, 4)
	R[j] = 0
			}
# Pacific Region 
	pc = ann4_med %*% S_pc
	reg_t  = lm(pc ~ tt)
	b_t = reg_t$coefficients[2]
	r_t = reg_t$residuals
	OMEGA0=omega0(r_t)
	VF = eta1 * b_t^2 / OMEGA0
	Table4_med[(loc_n+1),1] = round(b_t,4)
	Table4_med[(loc_n+1),2] = round( (b_t - 6.482*abs(b_t/chol(VF))), 4)
	Table4_med[(loc_n+1),3] = round( (b_t + 6.482*abs(b_t/chol(VF))), 4)
	Table4_med[(loc_n+1),4] = round(VF, 4)
	if (b_t > 0 & VF > 41.53) {Region[1,2] = " + "}
	if (b_t < 0 & VF < 41.53) {Region[1,2] = "(0)"}
	if (b_t < 0 & VF > 41.53) {Region[1,2] = " - "}
# SouthEast Region 
	se = ann4_med %*% S_se
	reg_t  = lm(se ~ tt)
	b_t = reg_t$coefficients[2]
	r_t = reg_t$residuals
	OMEGA0=omega0(r_t)
	VF = eta1 * b_t^2 / OMEGA0
	Table4_med[(loc_n+2),1] = round(b_t,4)
	Table4_med[(loc_n+2),2] = round( (b_t - 6.482*abs(b_t/chol(VF))), 4)
	Table4_med[(loc_n+2),3] = round( (b_t + 6.482*abs(b_t/chol(VF))), 4)
	Table4_med[(loc_n+2),4] = round(VF, 4)
	if (b_t > 0 & VF > 41.53) {Region[2,2] = " + "}
	if (b_t < 0 & VF < 41.53) {Region[2,2] = "(0)"}
	if (b_t < 0 & VF > 41.53) {Region[2,2] = " - "}
# All 
	all = ann4_med %*% S_all
	reg_t  = lm(all ~ tt)
	b_t = reg_t$coefficients[2]
	r_t = reg_t$residuals
	OMEGA0=omega0(r_t)
	VF = eta1 * b_t^2 / OMEGA0
	Table4_med[(loc_n+3),1] = round(b_t,4)
	Table4_med[(loc_n+3),2] = round( (b_t - 6.482*abs(b_t/chol(VF))), 4)
	Table4_med[(loc_n+3),3] = round( (b_t + 6.482*abs(b_t/chol(VF))), 4)
	Table4_med[(loc_n+3),4] = round(VF, 4)
	if (b_t > 0 & VF > 41.53) {Region[3,2] = " + "}
	if (b_t < 0 & VF < 41.53) {Region[3,2] = "(0)"}
	if (b_t < 0 & VF > 41.53) {Region[3,2] = "  -"}
#
cat("TABLE 4 (MED) COMPLETE", "\n")

#	VAR
reg_t  = lm(ann5_var ~ tt)					# regress annual stat (all locations) on time trend
b_t = reg_t$coefficients[2,]					# get slope coeffs
r_t = reg_t$residuals[,]					# Get residuals  
OMEGA1=omega1(r_t)
R=matrix(0, nrow=1, ncol=loc_n)
for (j in 1:loc_n)	{					# j-loop for locations
	R[j] = 1
	VF = eta1 * (R %*% b_t) %*% (solve( R %*% OMEGA1 %*% t(R) ) ) %*% (R %*% b_t)
	Table5_var[j,1] = round(b_t[j], 4)
	Table5_var[j,2] = round(b_t[j] - 6.482*abs( b_t[j]/chol(VF)), 4)
	Table5_var[j,3] = round(b_t[j] + 6.482*abs( b_t[j]/chol(VF)), 4)
	Table5_var[j,4] = round(VF, 4)
	R[j] = 0
			}
# Pacific Region 
	pc = ann5_var %*% S_pc
	reg_t  = lm(pc ~ tt)
	b_t = reg_t$coefficients[2]
	r_t = reg_t$residuals
	OMEGA0=omega0(r_t)
	VF = eta1 * b_t^2 / OMEGA0
	Table5_var[(loc_n+1),1] = round(b_t,4)
	Table5_var[(loc_n+1),2] = round( (b_t - 6.482*abs(b_t/chol(VF))), 4)
	Table5_var[(loc_n+1),3] = round( (b_t + 6.482*abs(b_t/chol(VF))), 4)
	Table5_var[(loc_n+1),4] = round(VF, 4)
	if (b_t > 0 & VF > 41.53) {Region[1,3] = " + "}
	if (b_t < 0 & VF < 41.53) {Region[1,3] = "(0)"}
	if (b_t < 0 & VF > 41.53) {Region[1,3] = " - "}
# SouthEast Region
	se = ann5_var %*% S_se
	reg_t  = lm(se ~ tt)
	b_t = reg_t$coefficients[2]
	r_t = reg_t$residuals
	OMEGA0=omega0(r_t)
	VF = eta1 * b_t^2 / OMEGA0
	Table5_var[(loc_n+2),1] = round(b_t,4)
	Table5_var[(loc_n+2),2] = round( (b_t - 6.482*abs(b_t/chol(VF))), 4)
	Table5_var[(loc_n+2),3] = round( (b_t + 6.482*abs(b_t/chol(VF))), 4)
	Table5_var[(loc_n+2),4] = round(VF, 4)
	if (b_t > 0 & VF > 41.53) {Region[2,3] = " + "}
	if (b_t < 0 & VF < 41.53) {Region[2,3] = "(0)"}
	if (b_t < 0 & VF > 41.53) {Region[2,3] = " - "}
# All 
	all = ann5_var %*% S_all
	reg_t  = lm(all ~ tt)
	b_t = reg_t$coefficients[2]
	r_t = reg_t$residuals
	OMEGA0=omega0(r_t)
	VF = eta1 * b_t^2 / OMEGA0
	Table5_var[(loc_n+3),1] = round(b_t,4)
	Table5_var[(loc_n+3),2] = round( (b_t - 6.482*abs(b_t/chol(VF))), 4)
	Table5_var[(loc_n+3),3] = round( (b_t + 6.482*abs(b_t/chol(VF))), 4)
	Table5_var[(loc_n+3),4] = round(VF, 4)
	if (b_t > 0 & VF > 41.53) {Region[3,3] = " + "}
	if (b_t < 0 & VF < 41.53) {Region[3,3] = "(0)"}
	if (b_t < 0 & VF > 41.53) {Region[3,3] = "  -"}
#
cat("TABLE 5 (VAR) COMPLETE", "\n")

#	MAX
reg_t  = lm(ann6_max ~ tt)					# regress annual stat (all locations) on time trend
b_t = reg_t$coefficients[2,]					# get slope coeffs
r_t = reg_t$residuals[,]					# Get residuals  
OMEGA1=omega1(r_t)
R=matrix(0, nrow=1, ncol=loc_n)
for (j in 1:loc_n)	{					# j-loop for locations
	R[j] = 1
	VF = eta1 * (R %*% b_t) %*% (solve( R %*% OMEGA1 %*% t(R) ) ) %*% (R %*% b_t)
	Table6_max[j,1] = round(b_t[j], 4)
	Table6_max[j,2] = round(b_t[j] - 6.482*abs( b_t[j]/chol(VF)), 4)
	Table6_max[j,3] = round(b_t[j] + 6.482*abs( b_t[j]/chol(VF)), 4)
	Table6_max[j,4] = round(VF, 4)
	R[j] = 0
			}
# Pacific Region 
	pc = ann6_max %*% S_pc
	reg_t  = lm(pc ~ tt)
	b_t = reg_t$coefficients[2]
	r_t = reg_t$residuals
	OMEGA0=omega0(r_t)
	VF = eta1 * b_t^2 / OMEGA0
	Table6_max[(loc_n+1),1] = round(b_t,4)
	Table6_max[(loc_n+1),2] = round( (b_t - 6.482*abs(b_t/chol(VF))), 4)
	Table6_max[(loc_n+1),3] = round( (b_t + 6.482*abs(b_t/chol(VF))), 4)
	Table6_max[(loc_n+1),4] = round(VF, 4)
	if (b_t > 0 & VF > 41.53) {Region[1,4] = " + "}
	if (b_t < 0 & VF < 41.53) {Region[1,4] = "(0)"}
	if (b_t < 0 & VF > 41.53) {Region[1,4] = " - "}
# SouthEast Region 
	se = ann6_max %*% S_se
	reg_t  = lm(se ~ tt)
	b_t = reg_t$coefficients[2]
	r_t = reg_t$residuals
	OMEGA0=omega0(r_t)
	VF = eta1 * b_t^2 / OMEGA0
	Table6_max[(loc_n+2),1] = round(b_t,4)
	Table6_max[(loc_n+2),2] = round( (b_t - 6.482*abs(b_t/chol(VF))), 4)
	Table6_max[(loc_n+2),3] = round( (b_t + 6.482*abs(b_t/chol(VF))), 4)
	Table6_max[(loc_n+2),4] = round(VF, 4)
	if (b_t > 0 & VF > 41.53) {Region[2,4] = " + "}
	if (b_t < 0 & VF < 41.53) {Region[2,4] = "(0)"}
	if (b_t < 0 & VF > 41.53) {Region[2,4] = " - "}
# All 
	all = ann6_max %*% S_all
	reg_t  = lm(all ~ tt)
	b_t = reg_t$coefficients[2]
	r_t = reg_t$residuals
	OMEGA0=omega0(r_t)
	VF = eta1 * b_t^2 / OMEGA0
	Table6_max[(loc_n+3),1] = round(b_t,4)
	Table6_max[(loc_n+3),2] = round( (b_t - 6.482*abs(b_t/chol(VF))), 4)
	Table6_max[(loc_n+3),3] = round( (b_t + 6.482*abs(b_t/chol(VF))), 4)
	Table6_max[(loc_n+3),4] = round(VF, 4)
	if (b_t > 0 & VF > 41.53) {Region[3,4] = " + "}
	if (b_t < 0 & VF < 41.53) {Region[3,4] = "(0)"}
	if (b_t < 0 & VF > 41.53) {Region[3,4] = "  -"}
#
cat("TABLE 6 (MAX) COMPLETE", "\n")



#	DPR
if (seasonal==1) {ann7_dpr[1,3]=1}				# If fall: give Fresno 1 dpr in 1889 so routine doean't crash
reg_t  = lm(ann7_dpr ~ tt)					# regress annual stat (all locations) on time trend
b_t = reg_t$coefficients[2,]					# get slope coeffs
r_t = reg_t$residuals[,]					# Get residuals  
OMEGA1=omega1(r_t)
R=matrix(0, nrow=1, ncol=loc_n)
for (j in 1:loc_n)	{					# j-loop for locations
	R[j] = 1
	VF = eta1 * (R %*% b_t) %*% (solve( R %*% OMEGA1 %*% t(R) ) ) %*% (R %*% b_t)
	Table7_dpr[j,1] = round(b_t[j], 4)
	Table7_dpr[j,2] = round(b_t[j] - 6.482*abs( b_t[j]/chol(VF)), 4)
	Table7_dpr[j,3] = round(b_t[j] + 6.482*abs( b_t[j]/chol(VF)), 4)
	Table7_dpr[j,4] = round(VF, 4)
	R[j] = 0
			}
# Pacific Region
	pc = ann7_dpr %*% S_pc
	reg_t  = lm(pc ~ tt)
	b_t = reg_t$coefficients[2]
	r_t = reg_t$residuals
	OMEGA0=omega0(r_t)
	VF = eta1 * b_t^2 / OMEGA0
	Table7_dpr[(loc_n+1),1] = round(b_t,4)
	Table7_dpr[(loc_n+1),2] = round( (b_t - 6.482*abs(b_t/chol(VF))), 4)
	Table7_dpr[(loc_n+1),3] = round( (b_t + 6.482*abs(b_t/chol(VF))), 4)
	Table7_dpr[(loc_n+1),4] = round(VF, 4)
	if (b_t > 0 & VF > 41.53) {Region[1,5] = " + "}
	if (b_t < 0 & VF < 41.53) {Region[1,5] = "(0)"}
	if (b_t < 0 & VF > 41.53) {Region[1,5] = " - "}
# SouthEast Region 
	se = ann7_dpr %*% S_se
	reg_t  = lm(se ~ tt)
	b_t = reg_t$coefficients[2]
	r_t = reg_t$residuals
	OMEGA0=omega0(r_t)
	VF = eta1 * b_t^2 / OMEGA0
	Table7_dpr[(loc_n+2),1] = round(b_t,4)
	Table7_dpr[(loc_n+2),2] = round( (b_t - 6.482*abs(b_t/chol(VF))), 4)
	Table7_dpr[(loc_n+2),3] = round( (b_t + 6.482*abs(b_t/chol(VF))), 4)
	Table7_dpr[(loc_n+2),4] = round(VF, 4)
	if (b_t > 0 & VF > 41.53) {Region[2,5] = " + "}
	if (b_t < 0 & VF < 41.53) {Region[2,5] = "(0)"}
	if (b_t < 0 & VF > 41.53) {Region[2,5] = " - "}
# All 
	all = ann7_dpr %*% S_all
	reg_t  = lm(all ~ tt)
	b_t = reg_t$coefficients[2]
	r_t = reg_t$residuals
	OMEGA0=omega0(r_t)
	VF = eta1 * b_t^2 / OMEGA0
	Table7_dpr[(loc_n+3),1] = round(b_t,4)
	Table7_dpr[(loc_n+3),2] = round( (b_t - 6.482*abs(b_t/chol(VF))), 4)
	Table7_dpr[(loc_n+3),3] = round( (b_t + 6.482*abs(b_t/chol(VF))), 4)
	Table7_dpr[(loc_n+3),4] = round(VF, 4)
	if (b_t > 0 & VF > 41.53) {Region[3,5] = " + "}
	if (b_t < 0 & VF < 41.53) {Region[3,5] = "(0)"}
	if (b_t < 0 & VF > 41.53) {Region[3,5] = "  -"}
#
cat("TABLE 7 (DPR) COMPLETE", "\n")



#	EXC 1 DAY, 99
reg_t  = lm(ann_1d99 ~ tt)					# regress annual stat (all locations) on time trend
b_t = reg_t$coefficients[2,]					# get slope coeffs
r_t = reg_t$residuals[,]					# Get residuals  
OMEGA1=omega1(r_t)
R=matrix(0, nrow=1, ncol=loc_n)
for (j in 1:loc_n)	{					# j-loop for locations
	R[j] = 1
	VF = eta1 * (R %*% b_t) %*% (solve( R %*% OMEGA1 %*% t(R) ) ) %*% (R %*% b_t)
	Table8_1day_99[j,1] = round(b_t[j], 4)
	Table8_1day_99[j,2] = round(b_t[j] - 6.482*abs( b_t[j]/chol(VF)), 4)
	Table8_1day_99[j,3] = round(b_t[j] + 6.482*abs( b_t[j]/chol(VF)), 4)
	Table8_1day_99[j,4] = round(VF, 4)
	R[j] = 0
			}
# Pacific Region 
	pc = ann_1d99 %*% S_pc
	reg_t  = lm(pc ~ tt)
	b_t = reg_t$coefficients[2]
	r_t = reg_t$residuals
	OMEGA0=omega0(r_t)
	VF = eta1 * b_t^2 / OMEGA0
	Table8_1day_99[(loc_n+1),1] = round(b_t,4)
	Table8_1day_99[(loc_n+1),2] = round( (b_t - 6.482*abs(b_t/chol(VF))), 4)
	Table8_1day_99[(loc_n+1),3] = round( (b_t + 6.482*abs(b_t/chol(VF))), 4)
	Table8_1day_99[(loc_n+1),4] = round(VF, 4)
	if (b_t > 0 & VF > 41.53) {Region[1,6] = " + "}
	if (b_t < 0 & VF < 41.53) {Region[1,6] = "(0)"}
	if (b_t < 0 & VF > 41.53) {Region[1,6] = " - "}
# SouthEast Region
	se = ann_1d99 %*% S_se
	reg_t  = lm(se ~ tt)
	b_t = reg_t$coefficients[2]
	r_t = reg_t$residuals
	OMEGA0=omega0(r_t)
	VF = eta1 * b_t^2 / OMEGA0
	Table8_1day_99[(loc_n+2),1] = round(b_t,4)
	Table8_1day_99[(loc_n+2),2] = round( (b_t - 6.482*abs(b_t/chol(VF))), 4)
	Table8_1day_99[(loc_n+2),3] = round( (b_t + 6.482*abs(b_t/chol(VF))), 4)
	Table8_1day_99[(loc_n+2),4] = round(VF, 4)
	if (b_t > 0 & VF > 41.53) {Region[2,6] = " + "}
	if (b_t < 0 & VF < 41.53) {Region[2,6] = "(0)"}
	if (b_t < 0 & VF > 41.53) {Region[2,6] = " - "}
# All 
	all = ann_1d99 %*% S_all
	reg_t  = lm(all ~ tt)
	b_t = reg_t$coefficients[2]
	r_t = reg_t$residuals
	OMEGA0=omega0(r_t)
	VF = eta1 * b_t^2 / OMEGA0
	Table8_1day_99[(loc_n+3),1] = round(b_t,4)
	Table8_1day_99[(loc_n+3),2] = round( (b_t - 6.482*abs(b_t/chol(VF))), 4)
	Table8_1day_99[(loc_n+3),3] = round( (b_t + 6.482*abs(b_t/chol(VF))), 4)
	Table8_1day_99[(loc_n+3),4] = round(VF, 4)
	if (b_t > 0 & VF > 41.53) {Region[3,6] = " + "}
	if (b_t < 0 & VF < 41.53) {Region[3,6] = "(0)"}
	if (b_t < 0 & VF > 41.53) {Region[3,6] = "  -"}

#
cat("TABLE 8 (1 DAY 99% EXC) COMPLETE", "\n")



#	EXC 2 DAY, 99
reg_t  = lm(ann_2d99 ~ tt)					# regress annual stat (all locations) on time trend
b_t = reg_t$coefficients[2,]					# get slope coeffs
r_t = reg_t$residuals[,]					# Get residuals  
OMEGA1=omega1(r_t)
R=matrix(0, nrow=1, ncol=loc_n)
for (j in 1:loc_n)	{					# j-loop for locations
	R[j] = 1
	VF = eta1 * (R %*% b_t) %*% (solve( R %*% OMEGA1 %*% t(R) ) ) %*% (R %*% b_t)
	Table9_2day_99[j,1] = round(b_t[j], 4)
	Table9_2day_99[j,2] = round(b_t[j] - 6.482*abs( b_t[j]/chol(VF)), 4)
	Table9_2day_99[j,3] = round(b_t[j] + 6.482*abs( b_t[j]/chol(VF)), 4)
	Table9_2day_99[j,4] = round(VF, 4)
	R[j] = 0
			}
# Pacific Region
	pc = ann_2d99 %*% S_pc
	reg_t  = lm(pc ~ tt)
	b_t = reg_t$coefficients[2]
	r_t = reg_t$residuals
	OMEGA0=omega0(r_t)
	VF = eta1 * b_t^2 / OMEGA0
	Table9_2day_99[(loc_n+1),1] = round(b_t,4)
	Table9_2day_99[(loc_n+1),2] = round( (b_t - 6.482*abs(b_t/chol(VF))), 4)
	Table9_2day_99[(loc_n+1),3] = round( (b_t + 6.482*abs(b_t/chol(VF))), 4)
	Table9_2day_99[(loc_n+1),4] = round(VF, 4)
	if (b_t > 0 & VF > 41.53) {Region[1,7] = " + "}
	if (b_t < 0 & VF < 41.53) {Region[1,7] = "(0)"}
	if (b_t < 0 & VF > 41.53) {Region[1,7] = " - "}
# SouthEast Region 
	se = ann_2d99 %*% S_se
	reg_t  = lm(se ~ tt)
	b_t = reg_t$coefficients[2]
	r_t = reg_t$residuals
	OMEGA0=omega0(r_t)
	VF = eta1 * b_t^2 / OMEGA0
	Table9_2day_99[(loc_n+2),1] = round(b_t,4)
	Table9_2day_99[(loc_n+2),2] = round( (b_t - 6.482*abs(b_t/chol(VF))), 4)
	Table9_2day_99[(loc_n+2),3] = round( (b_t + 6.482*abs(b_t/chol(VF))), 4)
	Table9_2day_99[(loc_n+2),4] = round(VF, 4)
	if (b_t > 0 & VF > 41.53) {Region[2,7] = " + "}
	if (b_t < 0 & VF < 41.53) {Region[2,7] = "(0)"}
	if (b_t < 0 & VF > 41.53) {Region[2,7] = " - "}
# All 
	all = ann_2d99 %*% S_all
	reg_t  = lm(all ~ tt)
	b_t = reg_t$coefficients[2]
	r_t = reg_t$residuals
	OMEGA0=omega0(r_t)
	VF = eta1 * b_t^2 / OMEGA0
	Table9_2day_99[(loc_n+3),1] = round(b_t,4)
	Table9_2day_99[(loc_n+3),2] = round( (b_t - 6.482*abs(b_t/chol(VF))), 4)
	Table9_2day_99[(loc_n+3),3] = round( (b_t + 6.482*abs(b_t/chol(VF))), 4)
	Table9_2day_99[(loc_n+3),4] = round(VF, 4)
	if (b_t > 0 & VF > 41.53) {Region[3,7] = " + "}
	if (b_t < 0 & VF < 41.53) {Region[3,7] = "(0)"}
	if (b_t < 0 & VF > 41.53) {Region[3,7] = "  -"}
#
cat("TABLE 9 (2 DAY 99% EXC) COMPLETE", "\n", "\n")
######################################################################################################


print(Region)


######################################################################################################
##	GENERATE OUTPUT FILES
s_code = c("Winter", "Spring", "Summer", "Fall")
sink("logfiles/tables.txt")
	cat("======================================================================","\n","\n")
	cat("SUMMARY TABLES - LONG TERM PRECIP FOR 10 US LOCATIONS","\n")
	cat("START YEAR: ", startyear, "; SEASONAL", (seasonal == 1), "; SEASON (IF APPLICABLE)", s_code[season.select],  "\n")
	cat("======================================================================","\n","\n")
	cat("Regional Trend Summaries:","\n")
	cat(" + pos sig, 0 pos insig, (0) neg insig, - neg sig; _w weighted","\n","\n")
	print(Region)
	cat("======================================================================","\n","\n")
	cat("Table S1: Summary Stats for 1-day Totals","\n")
	print(Table1_1day_total)
	cat("======================================================================","\n","\n")
	cat("Table S2: Summary Stats for 2-day Totals","\n")
	print(Table2_2day_total)
	cat("======================================================================","\n","\n")
	cat("Table S3: Trend Info for Annual Averages","\n")
	print(Table3_avg)
	cat("======================================================================","\n","\n")
	cat("Table S4: Trend Info for Annual Non-Zero Medians","\n")
	print(Table4_med)
	cat("======================================================================","\n","\n")
	cat("Table S5: Trend Info for Annual Variances","\n")
	print(Table5_var)
	cat("======================================================================","\n","\n")
	cat("Table S6: Trend Info for Annual Maximum","\n")
	print(Table6_max)
	cat("======================================================================","\n","\n")
	cat("Table S7: Trend Info for Annual Downpour Occurrences","\n")
	print(Table7_dpr)
	cat("======================================================================","\n","\n")
	cat("Table S8: Trend Info for Annual 1-day Exceedances of 1-day 99th percentile","\n")
	print(Table8_1day_99)
	cat("======================================================================","\n","\n")
	cat("Table S9: Trend Info for Annual 2-day Exceedances of 2-day 99th percentile","\n")
	print(Table9_2day_99)
	cat("======================================================================","\n","\n")
sink(file=NULL)


write.table(as.data.frame(Table1_1day_total), file="logfiles/Table1.csv", sep=",")
write.table(as.data.frame(Table2_2day_total), file="logfiles/Table2.csv", sep=",")


######################################################################################################
##	DRAW GRAPHS

pdf('graphics/Fig3.avg.pdf')
	yy = rep(year, loc_n)				# year series repeated for each location
	xy2 = stack(as.data.frame(ann3_avg) )
	xy3 = cbind(xy2,yy)
	print(	xyplot(xy3$values~xy3$yy|factor(xy3$ind), as.table=TRUE, 
		type=c("p","spline"), col="gray85", col.line="black", npoints=50, GCV=TRUE, xlab="Year", ylim=c(0,8), ylab="", main="Annual Average Precip (mm/day)",))
dev.off()

pdf('graphics/FigS1.med.pdf')
	xy2 = stack(as.data.frame(ann4_med) )
	xy3 = cbind(xy2,yy)
	print(xyplot(xy3$values~xy3$yy|factor(xy3$ind), as.table=TRUE, 
		type=c("p","spline"), col="gray85", col.line="black", npoints=100,  GCV=TRUE, xlab="Year", ylab="", ylim=c(0,12), main="Annual Median Precip (mm/non-zero day)"))
dev.off()

pdf('graphics/FigS2.var.pdf')
	xy2 = stack(as.data.frame(ann5_var) )
	xy3 = cbind(xy2,yy)
	print(xyplot(xy3$values~xy3$yy|factor(xy3$ind), as.table=TRUE, 
		type=c("p","spline"), col="gray85", col.line="black", npoints=30,  GCV=TRUE, ylab="",  xlab="Year", ylim=c(0,220), main="Annual Variance in Precip"))
dev.off()

pdf('graphics/Fig4.max.pdf')
	xy2 = stack(as.data.frame(ann6_max) )
	xy3 = cbind(xy2,yy)
	print(xyplot(xy3$values~xy3$yy|factor(xy3$ind), as.table=TRUE, 
		type=c("p","spline"), col="gray85", col.line="black", npoints=30,  GCV=TRUE, xlab="Year", ylab="", ylim=c(0,160), main="Annual Maximum (mm/day)"))
dev.off()

pdf('graphics/FigS3.dpr.pdf')
	xy2 = stack(as.data.frame(ann7_dpr) )
	xy3 = cbind(xy2,yy)
	print(xyplot(xy3$values~xy3$yy|factor(xy3$ind), as.table=TRUE, 
		type=c("p","spline"), col="gray85", col.line="black", GCV=TRUE, xlab="Year", ylab="", ylim=c(0,36), main="Annual Downpours (>1 inch/day)"))
dev.off()

pdf('graphics/FigS4.1d99.pdf')
	xy2 = stack(as.data.frame(ann_1d99) )
	xy3 = cbind(xy2,yy)
	print(xyplot(xy3$values~xy3$yy|factor(xy3$ind), as.table=TRUE, 
		type=c("p","spline"), col="gray85", col.line="black", xlab="Year", ylab="", ylim=c(0,10), main="Annual 1 day 99% Exceedances"))
dev.off()

pdf('graphics/Fig5.2d99.pdf')
	xy2 = stack(as.data.frame(ann_2d99) )
	xy3 = cbind(xy2,yy)
	print(xyplot(xy3$values~xy3$yy|factor(xy3$ind), as.table=TRUE, 
		type=c("p","spline"), col="gray85", col.line="black", xlab="Year", ylab="", ylim=c(0,16), main="Annual 2 day 99% Exceedances"))
dev.off()

#	DRAW RETURN EVENT GRAPHS
if(pentad.graph==1) 	{
y5=seq(startyear+4, endyear, 5)
c5 = c(1:length(y5))
p5pc = p5 %*% S_pc
p5se = p5 %*% S_se
pcb = lm(p5pc~c5)$coef[2]
seb = lm(p5se~c5)$coef[2]
print(summary(lm(p5pc~c5)))
print(summary(lm(p5se~c5)))
pdf('graphics/Figure6.pentad.pdf', width=10, height=6)	
plot(y5,p5pc, type="l", lwd=2, xlab="End Year of Pentad", ylab="Number of 1-in-5 Year Exceedances", col="blue", ylim=c(-1, 3), xaxt="n")
	abline(h=1)
	axis(side=1, at = y5)
par(new=TRUE)
	plot(y5,p5se, type="l", lwd=2, col="red", ylim=c(-1, 3), xlab="", ylab="", xaxt="n")
	text(1920, 3.0,"Pacific-blue ",col="blue", cex=1.2, pos=4)
	text(1957, 3.0, round(pcb,3), col="blue", pos=2)
	text(1920, 2.8,"SouthEast-red",col="red", cex=1.2, pos=4)
	text(1957, 2.8, round(seb,3), col="red", pos=2)
dev.off()

c50 = c(1943:2017)
p50pc = p50 %*% S_pc
p50se = p50 %*% S_se
pcb = lm(p50pc~c50)$coef[2]
seb = lm(p50se~c50)$coef[2]
print(summary(lm(p50pc~c50)))
print(summary(lm(p50se~c50)))
pdf('graphics/Figure6b.pentad50.pdf', width=10, height=6)	
plot(c50,p50pc, type="l", lwd=2, xlab="End Year of 50 yr Interval", ylab="Number of 1-in-50 Year Exceedances", col="blue", ylim=c(25, 40), xaxt="n")
	abline(h=30)
	axis(side=1, at = y5)
par(new=TRUE)
	plot(c50,p50se, type="l", lwd=2, col="red", ylim=c(25, 40), xlab="", ylab="", xaxt="n")
	text(1947, 40,"Pacific-blue ",col="blue", cex=1.2, pos=4)
	text(1967, 40, round(pcb,3), col="blue", pos=2)
	text(1947, 39,"SouthEast-red",col="red", cex=1.2, pos=4)
	text(1967, 39, round(seb,3), col="red", pos=2)
dev.off()
			}



######################################################################################################
##	ANALYSE LONG PROXY DATA SERIES & DRAW GRAPHS	
px = read.table("data/proxy.txt", header=TRUE)
pcp = px$PC
sep = px$SE
sep[sep == -99] = NA

#	DRAW TIME SERIES GRAPH
lma = 30
pcp_ma = filter(pcp, rep(1/lma,lma), sides=1)
	pcp_ma = pcp_ma - mean(pcp_ma, na.rm=TRUE) +1
sep_ma = filter(sep, rep(1/lma,lma), sides=1)
	sep_ma = sep_ma - mean(sep_ma, na.rm=TRUE) -1
ma=data.frame(pcp_ma,sep_ma)
pdf('graphics/Fig1.proxy.pdf', width=10, height=5)	
plot(px$Year,pcp_ma, xlab="Year", ylab="Precipitation Proxy", type="p", col="blue", ylim=c(-2.5,2.5))
	lines(sep_ma, col="red", type="p")
	abline(h=0)
	abline(v=1900, lty="dashed")
	text(100, 1.9,"Pacific",col="blue", cex=1.2)
	text(100,-0.5,"SouthEast",col="red", cex=1.2)
dev.off()

#	DRAW ACF GRAPHS
lag=c(1:200)
refp=.9^lag
refn=-.9^lag

z = acf(pcp, lag.max=200)
	zp=(z$acf)
z = acf(sep, lag.max=200, na.action=na.pass)
	zs=(z$acf)
pdf('graphics/Fig7.acf.pdf', width=8, height=6)
par(mfrow = c(1, 2))						# 2-Panel Figure
plot(zp[2:200], main="Pacific", type="p", ylab="Autocovariance Coefficients", xlab="Lag", ylim=c(-0.1,0.1), pch=16, col="gray")
	abline(h=0)
	par(new="TRUE")
	plot(refp[2:200], xlab="", ylab="", ylim=c(-0.1,0.1), type="l", lty="dashed", lwd=1.5, xaxt="n", yaxt="n")
	par(new="TRUE")
	plot(refn[2:200], xlab="", ylab="", ylim=c(-0.1,0.1), type="l", lty="dashed", lwd=1.5, xaxt="n", yaxt="n")
plot(zs[2:200], main="SouthEast", type="p", ylab="", xlab="Lag", ylim=c(-0.1,0.1), pch=16, col="gray")
	abline(h=0)
	par(new="TRUE")
	plot(refp[2:200], xlab="", ylab="", ylim=c(-0.1,0.1), type="l", lty="dashed", lwd=1.5, xaxt="n", yaxt="n")
	par(new="TRUE")
	plot(refn[2:200], xlab="", ylab="", ylim=c(-0.1,0.1), type="l", lty="dashed", lwd=1.5, xaxt="n", yaxt="n")
dev.off()

#	DRAW AVG TREND GRAPH
trends = Table3_avg[1:loc_n,1]
pdf('graphics/Fig2.TrendMap.pdf', width=8, height=5)
colour = rep("orange3", times = loc_n)
colour[trends>-0.004] = "orange2"
colour[trends>-0.002] = "orange1"
colour[trends>-0.001] = "gray85"
colour[trends> 0.001] = "skyblue1"
colour[trends> 0.002] = "skyblue2"
colour[trends> 0.004] = "skyblue3"
colour[trends> 0.006] = "royalblue1"
colour[trends> 0.008] = "skyblue4"
map("usa")
par(new = "TRUE")
plot(lat~lon, xlim=c(-123.3,-68), ylim=c(25,50), col=colour, type="p", pch=15, cex=2, xlab="", ylab="", xaxt="n", yaxt="n")
legend("topright",	legend=c("<-0.004","-0.004--0.002","-0.002--0.001","-0.001-0.001","0.001-0.002","0.002-0.004","0.004-0.006","0.006-0.008",">0.008"),
			col=c("orange3","orange2","orange1","gray85","skyblue1","skyblue2","skyblue3","royalblue1","skyblue4"),
			pch=15, cex=0.9, bg="white")

dev.off()


if(fd==1)	{

#	COMPUTE FRACTIONAL DIFFERENCE PARAMETERS
cat("FRACTIONAL DIFFERENCE (d) PARAMETERS:", "\n")
dpcp = fdSperio(pcp)
	cat("PACIFIC PROXY:   ", round(dpcp$d,3), "; STD DEV =", round(dpcp$sd.as,3), "RATIO",round(dpcp$d/dpcp$sd.as,3), "\n")
dsep = fdSperio(sep[366:2018])
	cat("SOUTHEAST PROXY: ", round(dsep$d,3), "; STD DEV =", round(dsep$sd.as,3), "RATIO",round(dsep$d/dsep$sd.as,3), "\n")
d1 = daily %*% S_pc
d1[is.na(d1)]=0
dpd = fdSperio(d1)
	cat("PACIFIC DAILY:   ", round(dpd$d,3), "; STD DEV =", round(dpd$sd.as,3), "RATIO",round(dpd$d/dpd$sd.as,3), "\n")
d1 = daily %*% S_se
d1[is.na(d1)]=0
dsd = fdSperio(d1)
	cat("SOUTHEAST DAILY: ", round(dsd$d,3), "; STD DEV =", round(dsd$sd.as,3), "RATIO",round(dsd$d/dsd$sd.as,3), "\n")
dpa = fdSperio(pc_avg)
	cat("PACIFIC ANNUAL", round(dpa$d,3), "; STD DEV =", round(dpa$sd.as,3), "RATIO",round(dpa$d/dpa$sd.as,3), "\n")
dsa = fdSperio(se_avg)
	cat("SOUTHEAST ANNUAL", round(dsa$d,3), "; STD DEV =", round(dsa$sd.as,3), "RATIO",round(dsa$d/dsa$sd.as,3), "\n")

		}

cat( "\n", "End of Program", "\n")