##
##
##  Bonus-Malus Markov Chain Example 
##
##
##-------------------------------------------


#--- Define Matrix Power Function  (P^n)---
m.power <- function(P,n){
    if (n==0) {
        P1 = diag(1, dim(P)[1])
    } else if (n==1){
        P1=P
    } else if (n>1) {
        P1=P
        for (i in 1:(n-1)){ P1 = P1 %*% P }
    }
    return(P1)
}



#--- Bonus-Malus ---
lambda=2

  a0 = dpois(0,   lambda)
  a1 = dpois(1,   lambda)
  a2 = dpois(2,   lambda)
  a3 = 1-ppois(2, lambda)

  A = c(a0,a1,a2,a3,  a0,0,a1,a2+a3,  0,a0,0,a1+a2+a3,  0,0,a0,a1+a2+a3 )
  P = matrix(A, 4,4, byrow=T)
  P

  m.power(P,10)




lambda=1.5

  a0 = dpois(0,   lambda)
  a1 = dpois(1,   lambda)
  a2 = dpois(2,   lambda)
  a3 = 1-ppois(2, lambda)

  A = c(a0,a1,a2,a3,  a0,0,a1,a2+a3,  0,a0,0,a1+a2+a3,  0,0,a0,a1+a2+a3 )
  P = matrix(A, 4,4, byrow=T)
  P

  m.power(P,10)






lambda=1

  a0 = dpois(0,   lambda)
  a1 = dpois(1,   lambda)
  a2 = dpois(2,   lambda)
  a3 = 1-ppois(2, lambda)

  A = c(a0,a1,a2,a3,  a0,0,a1,a2+a3,  0,a0,0,a1+a2+a3,  0,0,a0,a1+a2+a3 )
  P = matrix(A, 4,4, byrow=T)
  P

  m.power(P,10)







lambda=.5

  a0 = dpois(0,   lambda)
  a1 = dpois(1,   lambda)
  a2 = dpois(2,   lambda)
  a3 = 1-ppois(2, lambda)

  A = c(a0,a1,a2,a3,  a0,0,a1,a2+a3,  0,a0,0,a1+a2+a3,  0,0,a0,a1+a2+a3 )
  P = matrix(A, 4,4, byrow=T)
  P

  m.power(P,10)




lambda=.3

  a0 = dpois(0,   lambda)
  a1 = dpois(1,   lambda)
  a2 = dpois(2,   lambda)
  a3 = 1-ppois(2, lambda)

  A = c(a0,a1,a2,a3,  a0,0,a1,a2+a3,  0,a0,0,a1+a2+a3,  0,0,a0,a1+a2+a3 )
  P = matrix(A, 4,4, byrow=T)
  P

  m.power(P,10)









lambda=2

          [,1]      [,2]      [,3]      [,4]
[1,] 0.1353353 0.2706706 0.2706706 0.3233236
[2,] 0.1353353 0.0000000 0.2706706 0.5939942
[3,] 0.0000000 0.1353353 0.0000000 0.8646647
[4,] 0.0000000 0.0000000 0.1353353 0.8646647
>   m.power(P,10)
            [,1]       [,2]      [,3]      [,4]
[1,] 0.002692326 0.01718469 0.1215737 0.8585493
[2,] 0.002690677 0.01718294 0.1215724 0.8585539
[3,] 0.002689605 0.01718111 0.1215719 0.8585573
[4,] 0.002688979 0.01718066 0.1215714 0.8585590




lambda=1.5

          [,1]      [,2]      [,3]      [,4]
[1,] 0.2231302 0.3346952 0.2510214 0.1911532
[2,] 0.2231302 0.0000000 0.3346952 0.4421746
[3,] 0.0000000 0.2231302 0.0000000 0.7768698
[4,] 0.0000000 0.0000000 0.2231302 0.7768698
>   m.power(P,10)
           [,1]       [,2]      [,3]      [,4]
[1,] 0.01336314 0.04624358 0.1869593 0.7534340
[2,] 0.01331030 0.04620823 0.1869416 0.7535398
[3,] 0.01327269 0.04615059 0.1869487 0.7536280
[4,] 0.01324443 0.04614125 0.1869334 0.7536809




lambda=1

          [,1]      [,2]      [,3]      [,4]
[1,] 0.3678794 0.3678794 0.1839397 0.0803014
[2,] 0.3678794 0.0000000 0.3678794 0.2642411
[3,] 0.0000000 0.3678794 0.0000000 0.6321206
[4,] 0.0000000 0.0000000 0.3678794 0.6321206
m.power(P,10)
           [,1]      [,2]      [,3]      [,4]
[1,] 0.07256859 0.1223948 0.2590940 0.5459425
[2,] 0.07163789 0.1221962 0.2590628 0.5471031
[3,] 0.07090582 0.1212031 0.2596587 0.5482324
[4,] 0.07011132 0.1212655 0.2594601 0.5491631



lambda=.5

          [,1]      [,2]       [,3]       [,4]
[1,] 0.6065307 0.3032653 0.07581633 0.01438768
[2,] 0.6065307 0.0000000 0.30326533 0.09020401
[3,] 0.0000000 0.6065307 0.00000000 0.39346934
[4,] 0.0000000 0.0000000 0.60653066 0.39346934
m.power(P,10)
          [,1]      [,2]      [,3]      [,4]
[1,] 0.3743028 0.2396550 0.2082908 0.1777514
[2,] 0.3704601 0.2409249 0.2086774 0.1799376
[3,] 0.3670779 0.2369375 0.2134742 0.1825104
[4,] 0.3600241 0.2406091 0.2130136 0.1863531





lambda=.3
>   P
          [,1]      [,2]       [,3]        [,4]
[1,] 0.7408182 0.2222455 0.03333682 0.003599493
[2,] 0.7408182 0.0000000 0.22224547 0.036936313
[3,] 0.0000000 0.7408182 0.00000000 0.259181779
[4,] 0.0000000 0.0000000 0.74081822 0.259181779
>
>   m.power(P,10)
          [,1]      [,2]      [,3]       [,4]
[1,] 0.6246146 0.2176636 0.1067727 0.05094915
[2,] 0.6225156 0.2188820 0.1069689 0.05163351
[3,] 0.6205753 0.2167331 0.1101775 0.05251411
[4,] 0.6139414 0.2214267 0.1100187 0.05461316




#-------------------------------------
lambda=1

  a0 = dpois(0,   lambda)
  a1 = dpois(1,   lambda)
  a2 = dpois(2,   lambda)
  a3 = 1-ppois(2, lambda)

  A = c(a0,a1,a2,a3,  a0,0,a1,a2+a3,  0,a0,0,a1+a2+a3,  0,0,a0,a1+a2+a3 )
  P = matrix(A, 4,4, byrow=T)


D1 <- matrix(0, 16, 4)
for (i in 1:16){
    if (i==1) {
        D1[i,] <- c(0,0,0,1)
    } else {
        D1[i,] <- c(0,0,0,1) %*%  m.power(P,i-1)
    }
    plot(1:4, D1[i,], type="h", ylim=c(0,1), main=paste(i-1), lwd=12)
    readline("hit enter")
}

print(D1)

sum(  D1[16,] * c(200, 300, 400, 600) )









#-------------------------------------
#--- Coin toss Problem ---

A = c( 0,.5,.5, 0, 0, 0, 0, 0, 0,
       0, 0, 0,.5,.5, 0, 0, 0, 0,
       0, 0, 0, 0, 0,.5,.5, 0, 0,
       0, 0, 0,.5,.5, 0, 0, 0, 0,
       0, 0, 0, 0, 0,.5,.5, 0, 0,
       0, 0, 0, 0, 0,.5, 0,.5, 0,
       0, 0, 0, 0,.5, 0, 0, 0,.5,
       0, 0, 0, 0, 0, 0, 0, 1, 0,
       0, 0, 0, 0, 0, 0, 0, 0, 1  )

P = matrix(A, 9,9,byrow=T)


m.power(P,20)