# HW 3 Q 2 JLCY, 19 Dec 2013

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

# 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

# This code is to simulate from an ARMA fitted to the data for monthly
# streamflow time series.

# get the data into one long array...  scale the flow data.. - i.e. remove
# the monthly mean and divide by monthly standard deviation

X = LFmonflow.std

X = array(t(X))  #standardized anomalies..

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

# Fit various ARMA models p = order of AR component q = order of MA
# component

best_p = 0
best_q = 0
best_modelAIC = Inf

# search for least AIC with different p & q from 1 to 4 to find best p &
# best q
for (p in 1:4) {
    for (q in 1:3) {
        zz = arima0(X, order = c(p, 0, q))  # ARMAR model
        modelAIC = zz$aic  # find AIC
        if (best_modelAIC > modelAIC) {
            best_p = p
            best_q = q
            best_modelAIC = modelAIC
        }
    }
}

p = best_p
q = best_q

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

# Simulate

# ARMA with p = best_p, q = best_q
zz = arima0(X, order = c(p, 0, q))

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 = arima.sim(n = nmons, list(ar = c(zz$coef[1:p]), ma = c(zz$coef[(p + 
        1):(p + q)])), sd = sqrt(zz$sigma2)) + mean(X)
    # xsim=arima.sim(n=nmons,list(order=c(p,d,q),ar=c(zz$coef[1:2]),ma=c(zz$coef[3])),
    # sd=sqrt(zz$sigma2))+mean(X)

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

    # back standardize
    simdismon = simdismon * fsdev + fmean

    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")

}


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

par(ask = TRUE)

# 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")

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

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

##### Things you have to add/develop to this ###### Annual flow PDF Drought and
##### Surplus Stats for the annual flow First compute the annual flow - i.e. sum
##### the monthly flows Select Median of the observed annual flow as the
##### threshold Annual flows above this is surplus and below is drought Sum the
##### excess from the median for surplus and deficit from the median for the
##### deficit.  you can also compute the average/max length of drought and
##### surplus

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

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

## Max of obs drought =  9.144 
## Max of sim drought =  14.93 , Standard deivation =  2.057 
##  
## Max of obs surplus =  10.83 
## Max of sim surplus =  29.53 , Standard deivation =  2.134 
##  
## 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 =  34 , Standard deivation =  5.495 
##  
## Average length of osb 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.6313 
## 

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