'Computes the matrix :math:`P` and scalar :math:`d` that represent
the value function
.. math::
V(x) = x\' P x + d
in the infinite horizon case.  Also computes the control matrix
:math:`F` from :math:`u = - Fx`
Returns
P : array_like(float)
P is part of the value function representation of
:math:`V(x) = x\'Px + d`
F : array_like(float)
F is the policy rule that determines the choice of control
in each period.
d : array_like(float)
d is part of the value function representation of
:math:`V(x) = x\'Px + d`'

def stationary_values(self):

    (Q, R, A, B, N, C) = (self.Q, self.R, self.A, self.B, self.N, self.C)
    (A0, B0) = ((np.sqrt(self.beta) * A), (np.sqrt(self.beta) * B))
    P = solve_discrete_riccati(A0, B0, R, Q, N)
    S1 = (Q + (self.beta * dot(B.T, dot(P, B))))
    S2 = ((self.beta * dot(B.T, dot(P, A))) + N)
    F = solve(S1, S2)
    d = ((self.beta * np.trace(dot(P, dot(C, C.T)))) / (1 - self.beta))
    (self.P, self.F, self.d) = (P, F, d)
    return (P, F, d)
