# HW 3 Q 3 JLCY, 19 Dec 2013

# Gamma Distribution

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

##### Things you have to add/develop to this###### Annual flow PDF commands to
##### plot Annual flow stats

### The matrices of statistics are written in separate files You can
### use/modify the boxplots.R to create boxplots of the stats.  Again feel
### free to modify/change these codes.

## boxplot the May pdfs - commands at the end of this file Or you can write
## out the May PDFs from simulations and use/modify the boxpdfs.R file

# This code is to model a Seasonal AR(1) or Thomas-Fiering Model for monthly
# streamflow time series.

# you should HAVE THE DATE IN THIS MATRIX WYmonflow before you run this
# code..

# the data is in WYmonflow - that contains number of years as rows and 12
# columns (Jan through Dec)

# clear variables
rm(list = ls())

Qnumber = 12

source("C:\\Users\\jfklee90\\Desktop\\CVEN 6833\\HW 3\\loaddata.r")

## Package `sm', version 2.2-5: type help(sm) for summary information

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

# Simulate say ARMA(2,1)

# get the co-efficients of the Thomas Fiering Model gamma distribution from
# p.53 of Bras book

coef1 = 1:12
coef2 = 1:12
gammaej = 1:12

for (j in 1:12) {
    j1 = j - 1

    if (j == 1) 
        j1 = 12
    if (j == 1) 
        coef1[j] = (cor(LFmonflow[2:nyrs, j], LFmonflow[1:nyrs1, j1]))
    if (j > 1) 
        coef1[j] = cor(LFmonflow[1:nyrs, j], LFmonflow[1:nyrs, j1])

    gammaej[j] = (skew(LFmonflow[, j]) - ((coef1[j]^3) * skew(LFmonflow[, j1])))/(1 - 
        coef1[j]^2)^1.5

    coef2[j] = sqrt((var(LFmonflow[, j])) * (1 - coef1[j] * coef1[j]))
    coef1[j] = coef1[j] * sqrt(var(LFmonflow[, j])/var(LFmonflow[, j1]))
}

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

for (k in 1:nsim) {

    # number of values to be generated nyrs is the number of years and 12 is the
    # number of months.
    nmons = nyrs * 12

    xsim = 1:nmons

    i = round(runif(1, 1, nyrs))
    xsim[1] = LFmonflow[i, 1]
    xprev = xsim[1]

    for (i in 2:nmons) {
        j = i%%12
        if (j == 0) 
            j = 12

        j1 = j - 1
        if (j == 1) 
            j1 = 12

        # p.53 of Bras book
        x1 = xprev - mean(LFmonflow[, j1])
        x2 = coef2[j] * ((2/gammaej[j]) * (1 + (gammaej[j] * rnorm(1, 0, 1)/6) - 
            (((gammaej[j])^2)/36))^3 - (2/(gammaej[j])))

        xsim[i] = mean(LFmonflow[, j]) + x1 * coef1[j] + x2

        xprev = xsim[i]
    }

    # put the array inot a matrix
    simdismon = matrix(xsim, nrow = 12)
    simdismon = t(simdismon)

    maysim = simdismon[, 5]
    annsim = apply(simdismon, 1, sum)

    simpdf[k, ] = sm.density(maysim, eval.points = xeval, display = "none")$estimate
    annsimpdf[k, ] = sm.density(annsim, eval.points = annxeval, display = "none")$estimate

    source("C:\\Users\\jfklee90\\Desktop\\CVEN 6833\\HW 3\\sim_basicstat.r")

    source("C:\\Users\\jfklee90\\Desktop\\CVEN 6833\\HW 3\\sim_droughtstat.r")

}

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

# find out statistics
source("C:\\Users\\jfklee90\\Desktop\\CVEN 6833\\HW 3\\droughtstat.r")

## Max of obs drought =  9.144 
## Max of sim drought =  18.57 , Standard deivation =  2.35 
##  
## Max of obs surplus =  10.83 
## Max of sim surplus =  33.68 , Standard deivation =  2.749 
##  
## Sum of obs drought =  162.6 
## Sum of obs surplus =  202.8 
##  
## Average length of obs drought in years =  1.96 
##  
## Min length of obs drought in years =  1 
##  
## Max length of obs drought in years =  5 
## Max length of sim drought in years =  33 , Standard deivation =  5.399 
##  
## Average length of obs surplus in years =  2.042 
##  
## Min length of obs surplus in years =  1 
##  
## Max length of obs surplus in years =  6 
## Max length of sim surplus in years =  8 , Standard deivation =  0.6045 
## 

par(ask = TRUE)

# plot boxplots
source("C:\\Users\\jfklee90\\Desktop\\CVEN 6833\\HW 3\\boxplots6.r")

# plot may flow
source("C:\\Users\\jfklee90\\Desktop\\CVEN 6833\\HW 3\\mayflow.r")

# plot annual flow
source("C:\\Users\\jfklee90\\Desktop\\CVEN 6833\\HW 3\\annflow.r")

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