Skip to content
Merged
Changes from all commits
Commits
File filter

Filter by extension

Filter by extension

Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
56 changes: 30 additions & 26 deletions src/quadsums.jl
Original file line number Diff line number Diff line change
Expand Up @@ -7,37 +7,39 @@ Functions to compute quadratic sums
=#

@doc doc"""
Computes the expected discounted quadratic sum
var_quadratic_sum(A, C, H, beta, x_0)

Computes the expected discounted quadratic sum.

```math
q(x_0) = \mathbb{E} \sum_{t=0}^{\infty} \beta^t x_t' H x_t
```


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

##### Arguments
- `A::Union{Float64, Matrix{Float64}}` The `n x n` matrix described above (scalar)
if `n = 1`
- `C::Union{Float64, Matrix{Float64}}` The `n x n` matrix described above (scalar)
if `n = 1`
- `H::Union{Float64, Matrix{Float64}}` The `n x n` matrix described above (scalar)
if `n = 1`
- `beta::Float64`: Discount factor in `(0, 1)`
- `x_0::Union{Float64, Vector{Float64}}` The initial condtion. A conformable
array (of length `n`) or a scalar if `n = 1`
# Arguments

##### Returns
- `A::Union{Float64, Matrix{Float64}}`: The `n x n` matrix described above
(scalar if `n = 1`).
- `C::Union{Float64, Matrix{Float64}}`: The `n x n` matrix described above
(scalar if `n = 1`).
- `H::Union{Float64, Matrix{Float64}}`: The `n x n` matrix described above
(scalar if `n = 1`).
- `beta::Float64`: Discount factor in `(0, 1)`.
- `x_0::Union{Float64, Vector{Float64}}`: The initial condition. A conformable
array (of length `n`) or a scalar if `n = 1`.

- `q0::Float64` : Represents the value ``q(x_0)``
# Returns

##### Notes
- `q0::Float64`: Represents the value ``q(x_0)``.

# Notes

The formula for computing ``q(x_0)`` is ``q(x_0) = x_0' Q x_0 + v`` where

- ``Q`` is the solution to ``Q = H + \beta A' Q A`` and
- ``v = \frac{trace(C' Q C) \beta}{1 - \beta}``
- ``v = \frac{trace(C' Q C) \beta}{1 - \beta}``.

"""
function var_quadratic_sum(A::ScalarOrArray, C::ScalarOrArray, H::ScalarOrArray,
Expand All @@ -59,27 +61,29 @@ function var_quadratic_sum(A::ScalarOrArray, C::ScalarOrArray, H::ScalarOrArray,
end

@doc doc"""
Computes the quadratic sum
m_quadratic_sum(A, B; max_it=50)

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 `solve_discrete_lyapunov` for
doubling algorithm. See the documentation of `solve_discrete_lyapunov` for
more information.

##### Arguments
# Arguments

- `A::Matrix{Float64}` : An `n x n` matrix as described above. We assume in order
for convergence that the eigenvalues of ``A`` have moduli bounded by unity
- `B::Matrix{Float64}` : An `n x n` matrix as described above. We assume in order
for convergence that the eigenvalues of ``B`` have moduli bounded by unity
- `max_it::Int(50)` : Maximum number of iterations
- `A::Matrix{Float64}`: An `n x n` matrix as described above. We assume in order
for convergence that the eigenvalues of ``A`` have moduli bounded by unity.
- `B::Matrix{Float64}`: An `n x n` matrix as described above. We assume in order
for convergence that the eigenvalues of ``B`` have moduli bounded by unity.
- `max_it::Int(50)`: Maximum number of iterations.

##### Returns
# Returns

- `gamma1::Matrix{Float64}` : Represents the value ``V``
- `gamma1::Matrix{Float64}`: Represents the value ``V``.

"""
function m_quadratic_sum(A::Matrix, B::Matrix; max_it=50)
Expand Down