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

## Locfit smoothing of the 0,1 time series in time to look at the temporal
## variability of probability of being in state 1. 
### Locfit + GLM

yy = as.data.frame(cbind(ydes,timelen))
fit = locfit(ydes ~ timelen, data = yy, family="binomial")
plot(fit)
points(timelen[ydes==1],rep(0.4,length(timelen[ydes == 1])))
abline(h=0.5)


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)

#########

#Simulating from current year for two years
#
# suppose in the current year, say, 1999, the process was in state 1

# now simulate an ensemble of 100 states for the following two years.

statesim = array(NULL)
nsim = 100
for(i in 1:nsim){
xx = runif(1,0,1)
if(xx <= Pi[2,1])x1=0 else x1=1

if(x1 == 0){
xx = runif(1,0,1)
if(xx <= Pi[1,1])x2=0 else x2=1}

if(x1 == 1){
xx = runif(1,0,1)
if(xx <= Pi[1,1])x2=0 else x2=1}
statesim=rbind(statesim,c(x1,x2))

}

statesim=statesim[2:(nsim+1),]  #the first row will be NAs

### statesim will be a matrix of nsim rows and two columns. The columns will be simulated
## states 0, 1
## Thus you have an ensemmble forecast of states for 2000 and 2001.
## Now conditionally simulate the flow magnitude thus obtaining a matrix of the same
## size (nsim rows, 2 columns) - i.e., the ensemble forecast for the two years.

#YOu can compute the 5th and 95th percentiles of the columns and the median for the two years.
## repeat for all the years