function  absorbing(Q, R)
%
% This is a MATLAB function that calculates important characteristics 
% of an absorbing Markov chain.
% 
% The full Markov chain transition matrix P is assumed to be (k+m) by (k+m).
% The first k states are absorbing states; the last m states are 
% transient states.
% The matrix P is divided into 4 submatrices: k by k, k by m, m by k and m by m.
% The upper left submatrix is a k by k identity matrix (1's on the diagonal; 
% 0's elsewhere).
% The upper right submatrix is a k by m matrix of all 0's.
% The lower left submatrix is an m by k matrix R.
% The lower right submatrix is an m by m matrix Q.
% We input the matrices Q and R as data to the function.
% The matrix Q gives the transition probabilities among the m transient states.
% The matrix R gives the one-step transition probabilities from the m 
% transient states to the k absorbing states.
%
% Given the absorbing Markov chain characterized by the two matrices Q and R,
% we calculate three matrices describing the behavior of the Markov chain.
%
% We first echo the input:
Q
R
%
%
% The first matrix we calculate is the fundamental matrix N.
% The entry N(i,j) gives the expected number of visits to transient 
% state j before absorption, 
% starting in transient state i.
% The fundamental matrix N depends only on the square sub-transition matrix Q.
% It is not difficult to see that the total expected number of visits is the 
% sum of the expected number of
% visits on the different steps, and that the expected number of visits 
% to j on the n^th step, 
% starting in i
% is just Q(i,j).  Thus in matrix form, we have
% N = I + Q + Q^2 + Q^3 + Q^4 + ...  
%
% We now obtain a relatively simple formula for the fundamental matrix N.
% To do so, we observe that,
% if we multiply the left and right side by (I - Q), 
% then we obtain the equation N*(I-Q) = I
% because there is cancellation on the right side (it telescopes).
% Thus N = (I-Q)^{-1} = inv(I-Q),
% where inv is the matrix inverse function.
%
%
s = size(Q);
n = s(1);
I = eye(n);
N = inv(I - Q)
%
% Now we calculate the expected number of transitions before absorption, 
% starting from each transient state.
% We let m(i) be this expected number of transitions before absorption, 
% starting in state i.
% Thus m is an m by 1 column vector.  Clearly, m is made up of 
% the row sums of N.
%
w = ones(n,1);
m = N*w
%
% Now we calculate the probability of absorption in each of the k 
% absorbing states, 
% starting from each transient state.
% Let B(i,j) be the probability of being absorbed in absorbing state j 
% starting in transient state i.
% Note that 
%   B(i,j) = R(i,j) + (Q*R)(i,j) + (Q^2*R)(i,j) + (Q^3*R)(i,j) + ...
% so that B = (I + Q + Q^2 + Q^3 + ...)*R = N*R.
%
B = N*R
%
% We have computed the three matrices N, m and B.