#  Simulate a 2 state MC with transition matrix

alpha = .8
beta = .9

P = rbind( c(alpha, 1 - alpha), c(1 - beta, beta) )
x0 = 0

 MC2.sim = function (n = 100, x0 = 0, alpha = .8, beta = .9 ) 
{
 P = rbind( c(alpha, 1 - alpha), c(1 - beta, beta) )
 x = c(x0, rep(0,n))
 for( i in 1:n){
   if(x[i]==0){ pr = 1 - alpha}
   else{ pr = beta }
   x[i+1] = rbinom( 1, size =1, pr = pr)
 }
x
}

# count transitions from state i to j ,   i,j=0,1

MC2.counts = function (x) 
{
# x is a realization or sample path of 2 state MC
n = length(x)
# drop first element of x
y1 = x[-1]
# drop second element of x
y2 = x[c(1:(n-1))]

# num 0 -> 0
ind0 = y1==0
ind1 = y2==0
n.00 = sum(ind0 & ind1)

# num 0 -> 1
ind1 = y2==1
n.01 = sum(ind0 & ind1)

# num 1 -> 0
ind0 = y1==1
ind1 = y2==0
n.10 = sum(ind0 & ind1)

# num 1 -> 1
ind1 = y2==1
n.11 = sum(ind0 & ind1)

n.trans = list( n.00 = n.00  , n.01 = n.01 , n.10 = n.10 , n.11 = n.11)
y = cbind( y1 , y2)

# used for my debugging
#list(n.trans = n.trans , y = y)

n.trans

}


# simulate a single sample path for time t = 0 . 1 , ... , n
# and obtain the summary statistics

x0 = 0
n = 200
x.path = MC2.sim( n = n , x0 = 0 , alpha = .8 , beta = .7)

x.stats = MC2.counts(x.path)
x.stats

# from these the student can now calculate the maximum likelihood estimates for this sample path


# the student will now have to code simulate M iid sample paths and the resulting estimates for each path
# from this the student can then approximate the sampling distribution of sqrt(n)*(theta.hat - theta) using the idea of the parametric bootstrap



# what do sample paths look like for such a Markov chain?

x0 = 0
n = 50
x.path = MC2.sim( n = n , x0 = 0 , alpha = .8 , beta = .7)

x.stats = MC2.counts(x.path)
x.stats
plot( 0:n , x.path , type  = "l")