# Source code for quantecon.quadsums

"""
This module provides functions to compute quadratic sums of the form described
in the docstrings.

"""

import numpy as np
from numpy import sqrt, dot
import scipy.linalg
from .matrix_eqn import solve_discrete_lyapunov

[docs]def var_quadratic_sum(A, C, H, beta, x0):
r"""
Computes the expected discounted quadratic sum

.. math::

q(x_0) = \mathbb{E} \Big[ \sum_{t=0}^{\infty} \beta^t x_t' H x_t \Big]

Here :math:{x_t} is the VAR process :math:x_{t+1} = A x_t + C w_t
with :math:{x_t} standard normal and :math:x_0 the initial condition.

Parameters
----------
A : array_like(float, ndim=2)
The matrix described above in description.  Should be n x n
C : array_like(float, ndim=2)
The matrix described above in description.  Should be n x n
H : array_like(float, ndim=2)
The matrix described above in description.  Should be n x n
beta: scalar(float)
Should take a value in (0, 1)
x_0: array_like(float, ndim=1)
The initial condtion. A conformable array (of length n, or with
n rows)

Returns
-------
q0: scalar(float)
Represents the value :math:q(x_0)

Remarks: The formula for computing :math:q(x_0) is
:math:q(x_0) = x_0' Q x_0 + v
where

* :math:Q is the solution to :math:Q = H + \beta A' Q A, and
* :math:v = \frac{trace(C' Q C) \beta}{(1 - \beta)}

"""
# == Make sure that A, C, H and x0 are array_like == #

A, C, H = list(map(np.atleast_2d, (A, C, H)))
x0 = np.atleast_1d(x0)
# == Start computations == #
Q = scipy.linalg.solve_discrete_lyapunov(sqrt(beta) * A.T, H)
cq = dot(dot(C.T, Q), C)
v = np.trace(cq) * beta / (1 - beta)
q0 = dot(dot(x0.T, Q), x0) + v

return q0

[docs]def m_quadratic_sum(A, B, max_it=50):
r"""
Computes the quadratic sum

.. math::

V = \sum_{j=0}^{\infty} A^j B A^{j'}

V is computed by solving the corresponding discrete lyapunov
equation using the doubling algorithm.  See the documentation of
util.solve_discrete_lyapunov for more information.

Parameters
----------
A : array_like(float, ndim=2)
An n x n matrix as described above.  We assume in order for
convergence that the eigenvalues of :math:A have moduli bounded by
unity
B : array_like(float, ndim=2)
An n x n matrix as described above.  We assume in order for
convergence that the eigenvalues of :math:A have moduli bounded by
unity
max_it : scalar(int), optional(default=50)
The maximum number of iterations

Returns
========
gamma1: array_like(float, ndim=2)
Represents the value :math:V

"""

gamma1 = solve_discrete_lyapunov(A, B, max_it)

return gamma1