# HW 3 Q 7 JLCY, 19 Dec 2013

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

# clear variables
rm(list = ls())

library(sm)

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

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

nsim = 250

library(robust)

## Warning: package 'robust' was built under R version 3.0.2

## Loading required package: fit.models

## Warning: package 'fit.models' was built under R version 3.0.2

## Loading required package: lattice

## Warning: package 'lattice' was built under R version 3.0.2

## Loading required package: MASS

## Warning: package 'MASS' was built under R version 3.0.2

## 
## Attaching package: 'MASS'
## 
## The following object is masked from 'package:sm':
## 
##     muscle
## 
## Loading required package: robustbase

## Warning: package 'robustbase' was built under R version 3.0.2

## 
## Attaching package: 'robustbase'
## 
## The following object is masked from 'package:sm':
## 
##     aircraft
## 
## Loading required package: rrcov

## Warning: package 'rrcov' was built under R version 3.0.2

## Loading required package: pcaPP
## Loading required package: mvtnorm
## Scalable Robust Estimators with High Breakdown Point (version 1.3-4)

library(SMPracticals)

## Warning: package 'SMPracticals' was built under R version 3.0.2

## Loading required package: ellipse

## Warning: package 'ellipse' was built under R version 3.0.2

## 
## Attaching package: 'SMPracticals'
## 
## The following object is masked from 'package:robustbase':
## 
##     salinity
## 
## The following object is masked from 'package:MASS':
## 
##     cement, forbes, leuk, shuttle
## 
## The following object is masked from 'package:lattice':
## 
##     barley

library(HiddenMarkov)

## Warning: package 'HiddenMarkov' was built under R version 3.0.2

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

######### GLM ##### X contains

# Read in data LFann <-
# read.table('http://civil.colorado.edu/~balajir/CVEN6833/R-sessions/session3/files-4HW3/LF_1906-2005.txt')
LFann = read.table("C:\\Users\\jfklee90\\Desktop\\CVEN 6833\\R session 3\\4HW3\\LF_1906-2005.txt")

# ENSOann <-
# read.table('http://civil.colorado.edu/~balajir/CVEN6833/R-sessions/session3/files-4HW3/nino3_1905-2007.txt')
ENSOann = read.table("C:\\Users\\jfklee90\\Desktop\\CVEN 6833\\R session 3\\4HW3\\nino3_1905-2007.txt")

# PDOann <-
# read.table('http://civil.colorado.edu/~balajir/CVEN6833/R-sessions/session3/files-4HW3/PDO_1900-2011.txt')
PDOann = read.table("C:\\Users\\jfklee90\\Desktop\\CVEN 6833\\R session 3\\4HW3\\PDO_1900-2011.txt")

# AMOann <-
# read.table('http://civil.colorado.edu/~balajir/CVEN6833/R-sessions/session3/files-4HW3/AMO_1856-2011.txt')
AMOann = read.table("C:\\Users\\jfklee90\\Desktop\\CVEN 6833\\R session 3\\4HW3\\AMO_1856-2011.txt")

# Compile data from matching years
years <- seq(1906, 2005, 1)
# LFdata <-
# as.data.frame(cbind(LFann[,1],LFann[,2]/(10^6),ENSOann[match(years,ENSOann[,1]),2],PDOann[match(years,PDOann[,1]),2],AMOann[match(years,AMOann[,1]),2]))
# names(LFdata) <- c('year','LF','ENSO','PDO','AMO')

ydes = LFann[, 2]/(10^6)
Annual = LFann[, 2]/(10^6)
obsflow = ydes
N = length(ydes)

yth = median(LFann[, 2]/(10^6))
ydes[ydes < yth] = 0
ydes[ydes >= yth] = 1

ydes1 = ydes[-1]
ydeslag = ydes[-N]

ENSOann1 = ENSOann[match(years, ENSOann[, 1]), 2]
PDOann1 = PDOann[match(years, PDOann[, 1]), 2]
AMOann1 = AMOann[match(years, AMOann[, 1]), 2]

ENSOann1 = ENSOann1[-1]
PDOann1 = PDOann1[-1]
AMOann1 = AMOann1[-1]

# X <-
# as.data.frame(cbind(ydeslag,ENSOann[match(years,ENSOann[,1]),2],PDOann[match(years,PDOann[,1]),2],AMOann[match(years,AMOann[,1]),2]))

X = as.data.frame(cbind(ydeslag, ENSOann1, PDOann1, AMOann1))

names(X) <- c("LF1", "ENSO", "PDO", "AMO")

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

library(MASS)
# X = cbind(ydeslag,ENSOann1,PDOann1,AMOann1)
X = cbind(ydeslag, PDOann1)  # only important variable

glmmarkov = glm(ydes1 ~ X, family = binomial())

sim = predict(glmmarkov)

Nl = length(sim)
statesim = c(1:Nl)

for (i in 1:Nl) {
    pk = sim[i]
    statesim[i] <- ifelse(runif(1) < pk, 1, 0)
}

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

pair00 = matrix(0, length(ydes), 2)
pair01 = matrix(0, length(ydes), 2)
pair10 = matrix(0, length(ydes), 2)
pair11 = matrix(0, length(ydes), 2)

j1 = 1
j2 = 1
j3 = 1
j4 = 1

for (i in 1:(length(ydes) - 1)) {
    i1 = i + 1
    if (ydes[i] == 0) {
        if (ydes[i1] == 0) {
            pair00[j1, 1] = obsflow[i]
            pair00[j1, 2] = obsflow[i1]
            j1 = j1 + 1
        } else {
            pair01[j2, 1] = obsflow[i]
            pair01[j2, 2] = obsflow[i1]
            j2 = j2 + 1
        }
    } else {
        if (ydes[i1] == 0) {
            pair10[j3, 1] = obsflow[i]
            pair10[j3, 2] = obsflow[i1]
            j3 = j3 + 1
        } else {
            pair11[j4, 1] = obsflow[i]
            pair11[j4, 2] = obsflow[i1]
            j4 = j4 + 1
        }
    }
}

# remove 0s

row0_00 = apply(pair00, 1, function(row) all(row != 0))
pair00 = pair00[row0_00, ]

row0_01 = apply(pair01, 1, function(row) all(row != 0))
pair01 = pair01[row0_01, ]

row0_10 = apply(pair10, 1, function(row) all(row != 0))
pair10 = pair10[row0_10, ]

row0_11 = apply(pair11, 1, function(row) all(row != 0))
pair11 = pair11[row0_11, ]

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

par(ask = TRUE)

annflow = Annual
nsims = nsim

# statistics
armean = c(1:nsim)
arvar = c(1:nsim)
arcor = c(1:nsim)
arskw = c(1:nsim)
armax = c(1:nsim)
armin = c(1:nsim)

annxeval = seq(min(Annual) - 0.25 * sd(Annual), max(Annual) + 0.25 * sd(Annual), 
    length = 100)
xeval = annxeval

source("C:\\Users\\jfklee90\\Desktop\\CVEN 6833\\HW 3\\lag1-knn.03.R")

# a = lag1knn(Annual, nsim)

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

annsimpdf = simpdf

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

# boxplots

obsmean = mean(Annual)
obsvar = var(Annual)
obsskew = skew(Annual)
obslag1 = cor(Annual[1:99], Annual[2:100])
obsmax = max(Annual)
obsmin = min(Annual)

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

## Boxplot of statistics

# plot May
par(mfrow = c(2, 3))

boxplot(armean, main = "Mean")
par(new = T)
plot(obsmean, col = "red")

boxplot(arvar, main = "Var")
par(new = T)
plot(obsvar, col = "red")

boxplot(arskw, main = "Skew")
par(new = T)
plot(obsskew, col = "red")

boxplot(arcor, main = "Lag-1")
par(new = T)
plot(obslag1, col = "red")

boxplot(armax, main = "Max")
par(new = T)
plot(obsmax, col = "red")

boxplot(armin, main = "Min")
par(new = T)
plot(obsmin, col = "red")

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

## statesim will be a vector of 0s and 1s Thus you have an ensemble states
## Now conditionally simulate the flow magnitude for each state pair You need
## to write commands to pull the historical data into 4 state transition
## pairs 0,0; 0,1; 1,0; 1,1. Then use the flows from the appropriate
## transition pairs, Consider the pairs as Xt1 and Xt, find the K-NN of
## current flow in Xt1 and resample one of the neighbor and select the
## corresponding value of Xt as the simulated flow.

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