diff --git a/src/quadsums.jl b/src/quadsums.jl
index 07a87579..e42dd08c 100644
--- a/src/quadsums.jl
+++ b/src/quadsums.jl
@@ -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,
@@ -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)