# Loading packages and clearing workspace
#########################################
# install.packages("foreign")
# install.packages("isotone")
library(foreign) # Write files to .dta (or other) formats
library(isotone) # Isotonic regression, of which the pool-adjacent-violators algorithm is an example

rm(list=ls(all=TRUE))


# Pulling data from 4_export.do
###############################
# importing data
starsdata <- read.table("Data/Temp/export_to_R.txt", header=TRUE, sep=",") 
    # This relative file path will work if the R project is open



# Inspecting data and assigning series object names for efficiency
##################################################################
summary(starsdata)
attach(starsdata)


# Generating triangular matrices of input functions, used later to generate cumulative residuals
################################################################################################
n <- 1000
trireal <- matrix(data=NA, nrow=n, ncol=n)
for (i in 0:n) {
trireal[,i] <- c(rep(0,i), rep(1,n-i)*real[i])
}
tripsid <- matrix(data=NA, nrow=n, ncol=n)
for (i in 0:n) {
tripsid[,i] <- c(rep(0,i), rep(1,n-i)*psid[i])
}
trivhlss <- matrix(data=NA, nrow=n, ncol=n)
for (i in 0:n) {
trivhlss[,i] <- c(rep(0,i), rep(1,n-i)*vhlss[i])
}
triqog <- matrix(data=NA, nrow=n, ncol=n)
for (i in 0:n) {
triqog[,i] <- c(rep(0,i), rep(1,n-i)*qog[i])
}
tristudent <- matrix(data=NA, nrow=n, ncol=n)
for (i in 0:n) {
tristudent [,i] <- c(rep(0,i), rep(1,n-i)*student[i])
}
tricauchy_05 <- matrix(data=NA, nrow=n, ncol=n)
for (i in 0:n) {
  tricauchy_05 [,i] <- c(rep(0,i), rep(1,n-i)*cauchy_05[i])
}
tricauchy_15 <- matrix(data=NA, nrow=n, ncol=n)
for (i in 0:n) {
tricauchy_15 [,i] <- c(rep(0,i), rep(1,n-i)*cauchy_15[i])
}
tricauchy_2 <- matrix(data=NA, nrow=n, ncol=n)
for (i in 0:n) {
tricauchy_2 [,i] <- c(rep(0,i), rep(1,n-i)*cauchy_2[i])
}
tricSmooth <- matrix(data=NA, nrow=n, ncol=n)
for (i in 0:n) {
tricSmooth [,i] <- c(rep(0,i), rep(1,n-i)*cSmooth[i])
}

starsdata2<- data.frame(seq(1:n),xxx,cauchy_05,cauchy_15,cauchy_2,real,student,vhlss,psid,qog,cSmooth)
k <- 11


# Generating ratios of output densities to input densities
##########################################################
clist <- c("cb", "topJ", "rba", "rbnz", "minn", "explore_cb", "dataDriven_cb", "noExpData_cb", "pub_cb", "noPub_cb", "dataCode_cb", "dataCode_topJ", "eye_cb", "noEye_cb", "eye_topJ", "control")
clist2 <- c("cauchy_05", "cauchy_15", "cauchy_2", "student", "real", "vhlss", "psid", "qog", "cSmooth") 
for (i in clist){
	for (j in clist2){
	k <- k+1
	c <- get(paste("obs_",i,sep=""))/get(j)
	assign(paste("R_obs_",i,"_",j,sep=""),c)
	starsdata2[,k] <- get(paste("R_obs_",i,"_",j,sep=""))
	names(starsdata2)[k] <- paste("R_obs_",i,"_",j,sep="")
	}
}


# Calculating cumulative residuals from best increasing fit of above ratios
############################################################################
for (i in clist){
	for (j in clist2){
  print(paste("R_obs_",i,"_",j,sep=""))
  a <- get(paste("R_obs_",i,"_",j,sep=""))
  pavareg <- gpava(xxx, a, weights = get(j), solver = weighted.mean, ties = "primary", p = NA)
  c <- with(pavareg,x) # The x are fitted value outputs of the gpava function
  assign(paste("yf_",i,"_",j,sep=""),c)
  k <- k+1
	starsdata2[,k] <- get(paste("yf_",i,"_",j,sep=""))
	names(starsdata2)[k] <- paste("yf_",i,"_",j,sep="")
  res <- a - c
  temp <- get(paste("tri",j,sep=""))
  #print(paste("tri",j,sep=""))
  sumres <- temp%*%res
  sumy <- temp%*%a # This tracks the cdf for the observed z-scores (outputs) 
  d <- sumres/sumy[n] # Rescaling the residuals to read as a difference in densities
  assign(paste("S_",i,"_",j,sep=""),d)
  k <- k+1
	starsdata2[,k] <- get(paste("S_",i,"_",j,sep=""))
	names(starsdata2)[k] <- paste("S_",i,"_",j,sep="")
  # plot(xxx,d,"l",xlab="t-values", ylab="Partial sum of residuals", ylim=c(-0.01,0.04), xlim=c(0,10))
	}
}


# Passing new series back to STATA for graphing scripts
#######################################################
write.dta(starsdata2, "Data/Temp/export_from_r.dta", version = 11)