'Forecast the geometric sums
.. math::
S_x := E \Big[ \sum_{j=0}^{\infty} \beta^j x_{t+j} | x_t \Big]
S_y := E \Big[ \sum_{j=0}^{\infty} \beta^j y_{t+j} | x_t \Big]
Parameters
beta : scalar(float)
Discount factor, in [0, 1)
beta : array_like(float)
The term x_t for conditioning
Returns
S_x : array_like(float)
Geometric sum as defined above
S_y : array_like(float)
Geometric sum as defined above'

def geometric_sums(self, beta, x_t):

    I = np.identity(self.n)
    S_x = solve((I - (beta * self.A)), x_t)
    S_y = self.G.dot(S_x)
    return (S_x, S_y)
