# lag-1 2 X 2 transition matrix.
# Similarly for 3 X 3

#suppose the annual flow time series is in the vector, annflow

  medann = median(annflow)
 n = length(annflow)

#       - flows above this will be the wet state (1) and below this will
#        be dry state (0)

         states = 1:n
         states[annflow < medann] = 0
         states[annflow > medann] = 1 

	x = states

z = rep(0,n)

 pwet=length(x[x==1])/n
 pdry = length(x[xx=0])/n

for(k1 in 2:n){
	if(x[(k1-1)]==0 & x[k1]==0) z[k1] = 1
	if(x[(k1-1)]==0 & x[k1]==1) z[k1] = 2
	if(x[(k1-1)]==1 & x[k1]==0) z[k1] = 3
	if(x[(k1-1)]==1 & x[k1]==1) z[k1] = 4
}

n0=length(x[x==0])
if(x[n] == 0)n0=n0-1

n1=length(x[x==1])
if(x[n] ==1)n1=n1-1

p.00 = sum(z==1)/(n0)
p.11 = sum(z==4)/n1

p.01=1-p.00
p.10=1-p.11


#transition probability matrix..
Pi = matrix(c(p.00,p.10,p.01,p.11),2,2)


# simulate n years of state values from the markov chain..
xsim=1:n

#generate the first year..
uu=runif(1,0,1)
if(uu <= pwet)xsim[1]=1 else xsim[1]=0

xprev=xsim[1]

  for(i in 2:n){
  if(xprev==1)irow=2
 if(xprev==0)irow=1
 uu=runif(1,0,1)
  if(uu <= Pi[irow,1])xsim[i]=0 else xsim[i]=1
 xprev=xsim[i]
 }