Partial Solutions

We agreed that the candidate for \(T\) following the renewal-type argument given in the notes should be:

\[T_1 = \inf\{ n \geq 1: X_{n-1} =s, X_n=t \}\]

and that we should start the chain at \(X_0 = t\). Thus \(T_n\) would be the first time, we reach the states \(s\) and \(t\), in order, for the \(n\)th time. Breaking up the Markov chains by such return times allows us to use a regular law of large numbers to argue that \[S_{T_n}/T_n = n/T_n \to (\mathbb{E}T_1)^{-1} = \pi_s p_{st}.\]

We can check this fact, by recycling some old code, found in Number 2.

We enter the the transition matrix \(P\) into R

P <- matrix(c(1/4,1/4, 1/2,1/4,1/4,1/2,1/8,1/4,5/8), nrow =3)
P <-t(P)
P

##       [,1] [,2]  [,3]
## [1,] 0.250 0.25 0.500
## [2,] 0.250 0.25 0.500
## [3,] 0.125 0.25 0.625

The left-eigenvector corresponding to the stationary distribution is easily computed and normalized

eigen(t(P))

## eigen() decomposition
## $values
## [1]  1.000000e+00  1.250000e-01 -8.336361e-17
## 
## $vectors
##           [,1]          [,2]          [,3]
## [1,] 0.2752409  7.071068e-01  7.071068e-01
## [2,] 0.3853373  5.768999e-16 -7.071068e-01
## [3,] 0.8807710 -7.071068e-01  1.893603e-16

z=eigen(t(P))$vector[,1]
stat <- z/sum(z)
stat

## [1] 0.1785714 0.2500000 0.5714286

When we compute \(P^n\) for large \(n\), we expect to see the stationary measure repeated as a row vector in the matrix, since we know from our convergence theory on Markov chains that \(p_{ij}(n) \to \pi(j).\)

Q = P
for (i in 1:1000){
  Q <- Q %*% P
}
Q

##           [,1] [,2]      [,3]
## [1,] 0.1785714 0.25 0.5714286
## [2,] 0.1785714 0.25 0.5714286
## [3,] 0.1785714 0.25 0.5714286

In order to simulate the Markov chain, we can adapt the code we had used before:

step <- function(i){
  q = P[i,]
  x=-1
  u = runif(1)
  j=0
  cumq = cumsum(q)
  while(x==-1){
    j<-j+1
    if(u <= cumq[j]){x <-j}
  }
  x
}

steps <- function(y){
  x = y[1]
  n = y[2]
  for (i in 1:n){
    x <- c(x, step(x[i]))
  }
  x
}

We start the chain at \(3\) and simulates until it returns, with a while loop.

time <- function(i){
  n =0
  x=i
  n <- n+1
  x <-step(x)
  while(x !=3){
    n <-n+1
    x <- step(x)
  }
  n
}

m=mean(replicate(1000, time(3)))
1/m

## [1] 0.5646527

stat[3]

## [1] 0.5714286

We will verify the claim for the case \(s=1\) and \(t=2\).

time2 <- function(i){
  n =0
  x=i
  n <- n+1
  x = c(x, step(x[length(x)]) )
  while(x[length(x)-1] !=1 | x[length(x)] !=2   ){
    n <-n+1
    x = c(x, step(x[length(x)] ))
  }
  n
}

z=mean(replicate(10000, time2(2)))
1/z

## [1] 0.04395257

stat[1] * P[1,2]

## [1] 0.04464286

The following code here, was used to debug and make sure the above function worked properly–orginally it did not, since I mixed up and versus or.

time3 <- function(i){
  n =0
  x=i
  n <- n+1
  x = c(x, step(x[length(x)]) )
  while(x[length(x)-1] !=1 | x[length(x)] !=2   ){
    n <-n+1
    x = c(x, step(x[length(x)] ))
  }
  x
}

time3(2)

##  [1] 2 1 3 2 3 3 3 2 2 3 3 3 3 1 2

Endnotes

Version: 17 November 2024

Rmd Source

Return times

Extending what we learned

Partial Solutions

Endnotes