# HW 3 Q 5 JLCY, 19 Dec 2013

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

# PCA

# 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

# get variance matrix..
zs.arp = var(LFmonflow.std)

# do an Eigen decomposition..
zsvd.arp = svd(zs.arp)

# Principal Components... pcs.arp=t(t(zsvd.arp$u) %*% t(LFmonflow.std))
pcs.arp = LFmonflow.std %*% zsvd.arp$u

# Eigen Values.. - fraction variance
lambdas.arp = (zsvd.arp$d/sum(zsvd.arp$d))

par(mfrow = c(1, 1))
plot(1:7, lambdas.arp[1:7], type = "l", xlab = "Modes", ylab = "Frac. Var. explained")
points(1:7, lambdas.arp[1:7], col = "red")

par(ask = TRUE)

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

## Points to evaluate the May pdf

for (k in 1:nsim) {

    # simulation matrix - 12 rows and nyrs = # of years as columns
    simdismon = matrix(0, ncol = nyrs, nrow = 12)

    for (i in 1:12) {

        # only first 3 PCs important
        if (i <= 3) {
            arbest = ar(pcs.arp[, i], order.max = 25)
            p = order(arbest$aic)[1]
            zz = arima(pcs.arp[, i], order = c(p, 0, 0), include.mean = TRUE, 
                method = "ML")
            xsim = arima.sim(n = nyrs, list(ar = c(zz$coef[1:p])), sd = sqrt(zz$sigma2))
        } else if (i > 3) {
            xsim = rnorm(nyrs, mean = mean(pcs.arp[, i]), sd = sd(pcs.arp[, 
                i]))
        }
        simdismon[i, ] = xsim
    }

    simdismon1 = t(simdismon) %*% t(zsvd.arp$u)
    simdismon = t(t(simdismon1) * fsdev + fmean)

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

}

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

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

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

# 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 =  18.87 , Standard deivation =  2.464 
##  
## Max of obs surplus =  10.83 
## Max of sim surplus =  29.01 , Standard deivation =  1.988 
##  
## 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 =  46 , Standard deivation =  6.707 
##  
## 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.627 
##