MichaelJMath
6/2/2018 - 6:36 PM

Markov Process Function

Generate a markov process

def get_eigenvals(A):
    "Get eigenvalues of a matrix as a list"
    evects = A.eigenvects()
    result = [evects[i][0] for i in range(len(evects))]
    return result

def get_eigenvects(A):
    "Get eigenvectors of a matrix as a list"
    evects = A.eigenvects()
    results = [evects[i][2][0] for i in range(len(evects))]
    return results

def markov_process(transition_matrix, initial_state, iterations):
    """Get the state probabilities given a markov transition matrix of
    probabilities and an initial condition.
    
    Parameters
    ----------
    transition_matrix: sp.Matrix
        mxm matrix of transition probabilities
    initial_state: list, sp.Matrix
        mx1 vector of initial state values
    iterations: int
        number of iterations/steps to run of the process
    
    Returns
    --------
    pd.DataFrame
        contains probabilities of each state (columns) at each step (rows)
    """
    initial_state = sp.Matrix(initial_state)
    evects = get_eigenvects(transition_matrix)
    evals = get_eigenvals(transition_matrix)
    lams = sp.diag(*evals)
    S = sp.Matrix([x.T for x in evects]).T
    c = S.inv() * initial_state
    steady_state = c[0] * evects[0]
    result_dict = {}
    for i in range(iterations+1):
        result_dict[i] = np.array(S*(lams**i)*c).squeeze().astype(np.float64)
    idx_labels = ["State_"+str(i+1) for i in range(len(initial_state))]
    result_df = pd.DataFrame(result_dict, index=idx_labels).T
    result_df.loc['Steady_State', :] = np.array(steady_state).squeeze().astype(np.float64)
    return result_df