library(robust)
library(SMPracticals)
library(HiddenMarkov)

# Get Seasonal NINO3 DJF(1900) through SON(2009) - 440 values
nino3ses=scan("http://iridl.ldeo.columbia.edu/expert/SOURCES/.Indices/.nino/.EXTENDED/.NINO3/T/%28Dec%201899%29%28Nov%202009%29RANGE/T/3/boxAverage/data.ch")

## two states - create states you need..
## ydens is the state time series..

timelen = seq(1900,2010,by=0.25)
n1=length(timelen)-1
timelen=timelen[1:n1]

ydes = nino3ses
ydes[ydes >= 0]=1
ydes[ydes < 0]=0

fit = MClik(ydes)
plot(fit$order,fit$AIC)		#to select the best order

## note that the orders go from 0 to 3, so the first entry is for order 0

# suppose the best order is order 1
## need to manipulate

border=1

if(border == 1)Pi=fit$two

# compute the transition probability matrix
 rsums = rowSums(Pi)
 Pi = Pi %*% diag(1/rsums)

# unconditional probabilities of being in one of the two states (0,1)
nprob = fit$one / sum(fit$one)

## simulate from the transition probability
zsim = mchain(NULL,Pi,nprob)
zsim = simulate(zsim,nsim=N)


# see the tpm from the simulation
x=zsim
estPi <- table(x$mc[-length(x$mc)], x$mc[-1])
rowtotal <- estPi %*% matrix(1, nrow=nrow(Pi), ncol=1)
estPi <- diag(as.vector(1/rowtotal)) %*% estPi
print(estPi)

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

# Read in data
LFann <- read.table("http://civil.colorado.edu/~balajir/CVEN6833/R-sessions/session3/files-4HW3/LF_1906-2005.txt")
ENSOann <- read.table("http://civil.colorado.edu/~balajir/CVEN6833/R-sessions/session3/files-4HW3/nino3_1905-2007.txt")
PDOann <- read.table("http://civil.colorado.edu/~balajir/CVEN6833/R-sessions/session3/files-4HW3/PDO_1900-2011.txt")
AMOann <- read.table("http://civil.colorado.edu/~balajir/CVEN6833/R-sessions/session3/files-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)
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)

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

sim = predict(glmmarkov)

Nl = length(sim)
statesim=1:Nl

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

## statesim will be a vector of 0s and 1s
## Thus you have an ensemmble 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.