# Adam McCurdy's HW 3
setwd("C:/Users/Billy/Google Drive/CU Boulder Courses/CVEN 6833 Advanced Data Analysis/Homework/Homework 3")
source("adam_hw3_lib.R")

#Read in Lee's Ferry Monthly Flow
WYmonflow=matrix(scan("LeesFerry-monflows-1906-2006.txt"),ncol=13,byrow=T)
WYmonflow=WYmonflow[,2:13]
WYmonflow=WYmonflow/10^6             # convert to MAF

### Scale the data
fmean = apply(WYmonflow,2,mean)     #monthly mean
fsdev = apply(WYmonflow,2,sd)       #monthly standard deviation
wyArray = t((t(WYmonflow) - fmean) / fsdev)   # standardized anomalies matrix
wyArray=array(t(wyArray))       #standardized anomalies vector (in chronological order)

### Create ARMA model
armaGrid <- armaModeler(data = wyArray, 3, 3) # create table of possible ARMA models
bestARMA <- armaGrid[armaGrid$AIC == min(armaGrid$AIC),] # choose the best ARMA model based on the minimum AIC
armaModel <- arima(wyArray, order = c(bestARMA$P, 0, bestARMA$Q)) # apply best ARMA parameters to arima() to create model

### Simulate and transform simulated results based on observed mean and standard deviation
nsim=250    # number of simulations
simResults <- armaSimulator(nsim, armaModel, wyArray) # simulate best ARMA model to create representative flows
simResults <- backTransform(WYmonflow, simResults) # transform simulated results to reflect the mean and sd of the obs data

### Calculate statistics on simulated and observed data on monthly and annual timesteps 
armaStatsFrame <- armaStatsFun(simResults) # create data frame based on the monthly simulated ARMA model results with simulation number, month, statistical value
origStatsFrame <- armaStatsFun(WYmonflow) # create data frame based on the original monthly data simulation number, month, statistical value
annualStatsFrame <- annualStats(simResults) # create data frame for annual statistics of simulated ARMA model 
annualOrig <- annualStats(WYmonflow) # create data frame for annual statistics of observed data

### Create PDFs/boxplots of PDFs for simulated and observed results 
flowPDF(WYmonflow, simResults, 5, T) # create PDF for May flows
flowPDF(WYmonflow, simResults, "Annual", T) # create PDF for Annual flows

statsPlot(armaStatsFrame, origStatsFrame) # create monthly boxplots for ARMA and observed data statistics
annualPlot(annualStatsFrame, annualOrig) # create annual boxplots for ARMA and observed data statistics